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ABSTRACT 



This thesis describes the development of a computer program that 
calculates the propagation loss for low frequencies in a shallow ocean 
given the depth of source and receiver, the sound speed profile of the 
water, the frequency of the source, and the impedance and sound speed in 
the bottom. The program does this, by computing the sum of normal modes 
for a specified set of boundary conditions. At the surface, perfect 
pressure release is assumed, and the boundary condition at the bottom is 
one of impedance mismatch. An effort was made to develop a Fast Field 
Program, which would use a FFT to predict propagation loss at a variety 
of ranges by solving for a discrete set of wave numbers, but development 
was not completed. 



3 



TABLE OF CONTENTS 



I. INTRODUCTtON 8 

II. THEORY 13 

A. DERIVATION OF THE FIELD INTEGRAL 14 

B. NORMAL MODE METHOD FOR PRESSURE DETERMINATION 19 

C. DIRECT EVALUATION OF THE FIELD INTEGRAL 24 

D. ABSORPTION AND BOUNDARY EFFECTS 27 

III. DESCRIPTION OF PROGRAMS 33 

A. NORMAL MODE PROGRAM 34 

B. FAST FIELD PROGRAM 37 

IV. TESTS RUN ON PROGRAMS 39 

A. RIGID BOTTOM WITH NO ABSORPTION 40 

B. RIGID BOTTOM WITH ABSORPTION 51 

C. IMPEDANCE BOTTOM 63 

V. CONCLUSIONS AND RECOMMENDATIONS 78 

APPENDIX A MATRIX METHOD FOR FINDING HORIZONT.AL WAVE NUMBER 81 

APPENDIX 8 COMPUTER PROGRAM EXACT 83 

APPENDIX C COMPUTER PROGRAM FFP 91 

LIST OF REFERENCES 96 

INITIAL DISTRIBUTION LIST 97 



4 



LIST OF TABLES 



4.1 EVALUATIOtSl OF ACCURACY OF SUBROUTINE 'EIGENS' 41 

4.2 RESULTS FROM 'EIGENS' WITH GRADIENT AND WITH DUCTS 43 

4.3 PRESSURE COMPARISON - FFP VERSUS EXACT 59 

4.4 EXACT - NUMBER OF INCREMENTS AND CORRECTION FACTOR 66 

4.5 EXACT WITH IMPEDANCE BOTTOM 68 

4.6 EXACT SOLUTIONS VARYING CRITICAL ANGLES 70 



5 



LIST OF FIGURES 



1.1 Methods of Solving for the Pressure Field 11 

1.2 Solutions Derived from the Field Integral 12 

3.1 Depth Function Divergence 36 

4.1 Modal Pressure Variation with Depth 45 

4.2 Relative Stimulation of Real Modes 46 

4.3 Pressure Error Due to Wave Number Error 48 

4.4 Modal Pressures with Small Range Change 49 

4.5 Propagation Loss - No Absorption 50 

4.6 Interference Patterns for Modes 1 and 3 52 

4.7 Propagation Loss with Variable Absorption 53 

4.8 Propagation Loss for Sources at Varying Depths 54 

4.9 FFP Wave Number Pressure Contributions - No Absorption 56 

4.10 FFP Integral x^^ith Absorption 57 

4.11 FFP Wave Number Contributions - Varying Number of Samples ... 60’ 

4.12 FFP Integral with Varying Absorption 61 

4.13 Modal Pressure Function with Impedance Bottom 64 

4.14 Phase Shift/Bounce Versus Critical Angle 69 

4.15 Loss/Bounce and Loss/Meter Versus Critical Angle 72 

4.16 Mode 1 Pressure for Critical Angle 0.30 and 0.10 73 

4.17 Prop Loss - Critical Angle 0.1 Radian, 0.3 Radian 74 

4.18 Prop Loss Critical Angle 0.5 Radian, 0.9 Radian 75 

4.19 Prop Loss for Iso-Speed and Positive Gradient 77 



6 



ACKNOWLEDGEMENTS 



I am deeply Indebted to Dr. Tom Gabrielson whose teaching stimulated 
my interest in this topic and whose encouragement and direction made this 
thesis a meaningful experience. His leadership and example, together with 
a professional credo that knowledge involves much more than a familiarity 
with a library of complex formulae, inspire all who know him as a scien- 
tist or as a teacher. Dr. Alan Coppens proved a welcome source of knowl- 
edge and assistance. His patience and support were essential to the com- 
pletion of this thesis. 

I also wish to express my gratitude to the officers of the Canadian 
Armed Forces who made it possible for me to partake of this unique educa- 
tional and professional experience. 

Not least, I thank my wife, Birgit, for her unselfish support and 
understanding. 



7 



I. INTRODUCTION 



A wide variety of methods are available for determining the pressure 
field from a point source in a stratified ocean of constant depth. 
Figure 1.1 shows that, by using various transforms, it is possible to 
derive expressions describing the pressure field in terms of the linear 
wave equation, the Helmholtz equation and, finally, the field integral. 
Methods associated with these expresssions vary, and each has relative 
advantages and disadvantages in terms of speed, accuracy, and ability to 
deal with complex situations. Unfortunately, no single method is supe- 
rior in all respects and there must be trade-offs. In some situations, 
accuracy might be sacrificed for speed, or the ability to deal with 
unusual boundaries may override the importance of speed or accuracy. 

Figure 1.2 shows the common methods of solving the field integral to 
determine the pressure. The Normal Mode Solution, the Method of Stee- 
pest Descent, The Multiple Scattering Method and Direct Numerical Integ- 
ration are "exact" methods. The Fast Field Program (FFP), a method of 
numerical integration, must be considered an approximation, not only 
because it approximates a Hankel function by an exponential, but also 
because of its reliance on a finite number of samples. Since Stationary 
Phase normally requires approximations, then so does its derivative, the 
Method of Images. The Wentzel-Kramers-Brillouin (WKB) Integral normally 
involves approximations so both Ray Theory and the WKB Mode Equation 
will seldom be exact. It is possible to find exact solutions using the 
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Normal Mode Technique, but in this paper, as will be seen in Chapter II, 
practical assumptions relegate the solutions to the approximate cate- 
gory. This thesis is limited to consideration of the Normal Mode Solu- 
tion, which is computed by the program EXACT, and a Fast Field Program, 
computer program FFP. 

Chapter II derives the field integral, then develops the Normal Mode 
Solution, and describes the Fast Field Program, FFP. The basic differ- 
ence between the methods involves the requirement of the Normal Mode 
Solution for an accurate estimation of the horizontal wave numbers for 
each of the modes. The FFP does not calculate the modal wave numbers 
but computes the pressure contributions for a large number of pre- 
defined horizontal wave numbers. An expression to describe losses and 
phase shift in terms of the impedance and sound speed ratios between the 
bottom and the water is also derived in Chapter II. 

Chapter III describes the programs EXACT and FFP. EXACT solves the 

field Integral by first finding good approximations of the mode wave 

numbers using a matrix eigenvalue technique and then refining them. FFP 
calculates the pressure field at a discrete number of horizontal wave 
numbers without regard to modal values. Both programs use the vector 
sum of the components of pressure from the constituent wave numbers to 
calculate the propagation loss between the source and the receiver, 
taking into account absorption in the water, losses to the bottom, and 
phase changes at the surface and bottom. 

Chapter IV describes the testing procedures and the results. Within 
this chapter the testing criteria are defined and the adaptive nature of 

the development process is discussed. Test results provided some very 
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instructive lessons which gave an insight into how normal modes can 
describe pressure variations in both the vertical and horizontal. 

Chapter V discusses the results and some of the weaknesses of the 
progr am. 
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II. THEORY 



In this discussion of normal modes, frequency will be very low so 
the absorptive losses due to the water will be negligible, except at 
long ranges. However, absorption in the bottom will often be very high, 
especially for the higher modes. The theory was used to derive formulae 
for the lossless case; then absorption was introduced as an imaginary 
component of the horiztonal wave number. Although absorption in the 
water is low for low frequencies, it proved expedient to input abnor- 
mally high absorption rates into the test programs. By doing this, it 
was possible to verify that the complex numbers were producing consis- 
tent values. 

The theoretical approach involves a series of transforms, some of 
which require approximations for simplification. These transforms 
successively take the wave equation from the time domain to the fre- 
quency domain, from three dimensions to two, from range dependent to 
wave number dependent, and finally to a function that describes pressure 
as a function of the depths of the receiver and transmitter, the hori- 
zontal wave number and the horizontal range. Normal mode theory 
enables one to consider the total pressure as the sum of pressures from 
many components. Each mode is stimulated by the transmitter, or sensed 
by the receiver, to a degree that depends on the positions of the 
transmitter and receiver in relation to the modal depth function curve. 
This curve is described in Chapter IV. 
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In order to keep the program simple, several assumptions must be 
made; some limit the practicality of this particular program and could 
be overcome in a more comprehensive routine, while others are of little 
consequence. 

The main assumptions, some of which will be discussed in detail 
later, are: 

1. source is cw and monochromatic; it emits a continuous wave 
at a constant frequency: 

2. source is a point radiating uniformly in all directions; 

3. medium is homogeneous in all respects in the horizontal; 

4 . bottom and surface are flat and parallel; 

5. surface is perfectly reflecting: 

6. bottom can be completely described by its impedance; 

7. all energy transmitted into the bottom is lost, and 

8. branch line integrals can be ignored. 



A. THE FIELD IfITEGRAL 

An omnidirectional point source operating at frequency f(t) is posi- 
tioned at a depth z in a water mass that is uniformly homogeneous in 
sound speed in. the horizontal plane and can be described in the vertical 
plane by the sound speed, c(z), at any depth. 

As long as the variations are linear, the pressure is given by the 
acoustic wave equation: 



V^p - 



c^z)- 



dP 



-f(t)6(s“So) 



2.1 



where: 

p - p(s,t) pressure as a function of position and time, 
f(t) » function describing the source amplitude. 



lA 



c(z) = depth dependent sound speed, and 
s = position in Cartesian coordinates. 

A Fourier transform from time to frequency, to, corresponds to con- 
sidering the equation for distinct frequencies and yields: 

P = -F(to)6(s-So) 2.2 



where: 

P = P(s,to) a function of position and frequency, 

to = frequency in radians/sec, and 

F(to) = transform of the driving function, f(t). 

Because the sound field is independent of the azimuthal angle, this 
transform can be formulated in the more manageable cylindrical co-ordi- 
nate system. In the new system the wave equation appears as: 



8^P , 1 3P 

3r^ r 9r 



3z^ 



+ k^P 



-F(to) 



6(z-Zo) 

2iTr 



dr 



2.3 



where: 

k = wave number, — , 
c(z)’ 

z = depth, and 

P = P(r,z,to), a function of horizontal range, depth and frequency. 

The Hankel transform allows one to compute the changes of pressure 
with a change of depth, independent of the range. It requires consider- 
ation of the vertical component of the wave number, k^: 
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where: 



i!2 

9z^ 



+ kzp = 



_F_ 

2tt 



6(z-Zo) • 



2.H 



p = p(5,rw) a function of horizontal wave number, horizontal range and 
frequency, 

kz = 

5 = horizontal wave number, and 

k = wave number, f , . 

c(z) 

This equation is then solved by variation of parameters. Let 



p = Wi(z)u(z) + W 2 (z)v(z) 2.5 



where: 



w.U) - 



f z 






v(z)6(z-Zo)dz 



W, 



2.5a 



uv 



W,(z) = ^ 



F_ 

2 IT 



u(z)6(z-Zo)dz 



W. 



2.5b 



uv 



Wyy “ ^ ~ ^ ^ (the Wronskian), 2.5c 

and both u and v are solutions of the homogeneous form of equation 2.4. 

Mow if v(z) satisfies the upper boundary condition and u(z) satisfies 
the lower boundary condition, then: 
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w 



1 



(_tr v(Zo) 






2ir Wuv 
0 



z > Zo 

Z ^ Z Q 



2.6a 



and 



«2 



F u(Zq) 
2it Wuv 
0 






z < Zo 
z > Zo 



2.6b 



This means that the transformed pressure can be written in terms of depth 
dependent u(z) and v(z) as: 



P = 



_F_ 

2it 



u(Zo)v(z) 

Wuv 



for z < Zo 



2.7a 



P 



p u(z)v(Zo) 
^ Wuv 



for z > Zo 



2.7b 



A convenient way of representing this is 




u(z>)v(z<) 



2.8 



where: 

z> represents the larger of z and Zo, and 
z< represents the smaller of z and Zo- 

To obtain the actual pressure one must next take the inverse Hankel 
transform of zero order: 
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p = 



pJo(5r)CdC 



2.9 



_F_ 

2n 



u(z>)v(z<) 



W. 






2.10 



uv 



Use is made of the fact that [Ref. 1 ] 



Jo(Cr) = 1/2 






2.1 1 



and 






2.12 



It is then possible to express the field integral (Equation 2.10) in 
terms of a Hankel function of the first kind [Ref. 2]: 



P 



00 



IL 

4it 

; - 00 



u(z>)v(z<) 



H(5^\cr)CdC 



2.13 



The Hankel function can be approximated by 




Ur) 



J_2 ^l(5r-./i» ( 

V Ti^r 



+ 



_1 

18 5r 



+ 



2.14 
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Restricting this .expansion to the first term makes the integral much 
more manageable but restricts the minimum horizontal range so that ^r 
>> 1 : 



P 



_F_ 

4tt 




f 00 



-00 



u(z>)v(z<) 






2.15 



This field integral represents the pressure, within the limits 
imposed by truncating the expansion of the Hankel function. This integ- 
ral forms the basis for the derivation of formulae describing pressure 
using two different numerical approaches. At this point, the approach to 
the normal mode solution, which will be used in the program EXACT, takes 
a different tack from the approach which will be used in the FFP solu- 
tion and program. 



B. NORMAL MODE METHOD FOR PRESSURE DETERMINATION 

In this paper the normal mode method refers to the process whereby 
approximations for the eigenvalue of a mode are successively refined by 
numerically integrating the square of a depth function across the depth. 
A correction term is developed from this integral and is used to refine 
the approximation until the correction is negligible. 

Let the guess at the normal mode horizontal wave number, be 5. 
For v(z) satisfies 



V" + k| V = 0 



2.16a 
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where; 



k 



2 ■ _ 
Z “ 






For 5n> ^nCz) satisfies 



"zn 



V = 0 



2. 1 6b 



where 



k|n = 

Multiplying 2.1 6a by Vj^ and 2.16b by v and subtracting yields 



vnv" - vvn" + vvn(k|-k|n) = 0 



2.17 



Then 



dz 



(vnv'-vvn') 



+ vvn(k^C^~k^+Cn^) 



0, or 



,_d_ 

dz 



(vnv’-vvn’) = vvn. 



2.18 



Integrating from depth 0 to depth h: 



(vnV 




(C^ - 5^) 



vvn dz 



2.19 



or 
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Vj^(h)v'(h) - Vn'(h)v(h) - Vn(0)v'(0) + Vn(0)v(0) = (5^"Cn) 



h 



vvj^dz. 2.20 

Jo 

The following argument holds for any set of permissible boundary 
conditions but, in order to clarify the discussion, the specific case of 
pressure release at z = 0 and rigid at z = h will be selected. The 
function v(z) is not a mode, but it can be forced to satisfy one of the 
boundary conditions without losing generality. If v(z) ia defined to 
satisfy the boundary condition at z = 0, then, for any permissible boun- 
dary condition, the last pair of terms on the left hand side of equation 
2.20 is identically zero. Therefore, 



vrj(h)v’(h) - vn'(h)v(h) = ( 



vvj^ dz. 



2.21 



For the simplest case considered in this paper, the bottom boundary con- 
dition is rigid so that Vj^'(h) = 0. Therefore 



vn(h)v’(h) = (c" 









vVfj dz 

. 0 



2.22 



and, since for 5 close to 5^^, 



AC" 



v(h)v'(h) 

/vMz 
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- 



v(h)v^ (h) 
fv^dz 



2.23 



Besides proving useful in determining successively better approximations 
of this same integral forms part of the expression for the contribu- 
tion of a mode, since it is proportional to the derivative of the Wron- 
skian with respect to wave number. 

The Cauchy integral enables one to write equation 2.13 as the sum of 
the residues of all the poles of the integrand. Each pole corresponds 
to a normal mode. 



P 



iL 

2tt 












Cn 



2.24 



where: 



W(Cn) = 0 

To derive an expression for (dW/d^) multiply both sides by 

Un(b) Un' 

^ = Vh) ^ ^ 



Un(h)v'(h) - u/i(h)v(h) = (C~Cn)U+Cn)S 



vMz 



2.25 



The left-hand side of this equation gives the Wronskian for u and v, 
and, as 5 > it should: 
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dC 



2.26 




W(5)-W(Cn) 



8(5+5n) 



'h 

io 



v^dz 



2.27 



Since 5 _ and by properly constructing the program, the constant g 
can be forced to be equal to 1 , 



dC 



^n 



fh 



25n 



v^dz 

. 0 



2.28 



and equation 2.24 can now be written: 



P 




n 



^n^n 

2^i^fv^dz 



^0 ^(5n*")5n' 



2.29 



Using the exponential for the Hankel function, as in eqn. 2.13 = 



iL y ,/_2_ ^n^n i^nf' 
4n ^ V /v^dz 



2.30 




n 




UnVne 



Un*" 



i tt/4 



/vMz 



2.31 
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In reality, this only represents the pressure in a system for which 
there are no boundary losses: i.e., for self-adjoint boundary conditions. 
If real losses, including scattering, were to be taken into account. 



rh 



the term 



JO 



vMz must be replaced by [Ref. 3] = 



•h 
. 0 



vMz 



v^(0)(^) - v^(h)(^) 
Cic dc _ 



2.32 



where: 



, dv/dz du/dz 

A = u = 

V u 

Ignoring boundary losses, equation 2.31 represents the total pressure 
at a point which is at depth z and horizontal distance r from a point cw 
source at depth Zq. This equation is the basis of the computer program 
EXACT, which is discussed in detail in Chapter III. 

C. DIRECT EVALUATION OF THE FIELD INTEGRAL 

The second method of determining pressure at a point requires calcu- 
lating the field integral itself and then numerically integrating. Equa- 
tion 2.15 represents the pressure field in terms of the Wronskian and 
the source and receiver depth functions u(z) and v(z). The integral is 
solved by replacing it with a summation across an interval of wave 
numbers. Pressure at a point is found to be the sum of the pressure 
contributions from each horizontal wave number: 
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p 



’^mA£ A5 



2.33 



_F_ 

4tt 




-itt/4 

e 






“£n 



u(z>)v(z<) 



imA£r 

e 



u(z>), v(z<), and are calculated from the known boundary conditions 

at the surface and at the bottom. A£ must be small enough to ensure 
accuracy. 

It is impractical to pursue this theme without introducing absorp- 
tion. Evaluation of equation 2.3^ when mA£ = produces = 0 in the 
denominator, leading to an infinite solution. The magnitude of then 
depends on how near mA£ comes to 5^. Introduction of absorption as an 
imaginary component of horizontal wave number displaces the poles from 
the real axis and ensures that the Wronskian does not go to zero along 
the path of integration. This has the effect of dispersing the energy 
out of the theoretically limitless pressure function that results when 
the Wronskian goes to zero at the exact wave number corresponding to the 
normal modes. 

The terms within the summation sign of equation 2.33 have the same 
form as the discrete Fourier transform, normally considered for time, t, 
and frequency, oj: 

N-1 

G(nA(j) e^2imm/N ^ g(mAt). 2.3^ 

n=0 



By simply renaming the variables, the transform becomes: 
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X. G(nA5) = g(mAr). 2.35 

n=0 

This is the basis for the technique known as the Fast Field Program 
[Ref.^]. By comparing equation 2.3^ to equation 2.35, the following can 
be seen to be true: 



G(nA^) 



uv 

W 



5 






The pressure at range R can then be written as 



2.36 



^R “ ^mAr 



_F^ 

^7T 



/ 



TTtnAr 



-iir/4 

e 



A5[FFT[G]]. 



2.37 



The capability of this Fast Field Program will be limited by the number 
of range (and wave number) increments. The wave number sampling incre- 
ment, A5, the range resolution, Ar, and the number of points in the FFT 
can be seen to be related by: 



2Trnm 

N 



Cr = (nA^)(mAr) 



2.38 



or 



2n 

N 



A^Ar, 



2.39 
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which results from comparing equations 2.3^ and 2.35. The number of 
sampling points, N, places practical limits on either the maximum range 
or the range resolution. 



D. ABSORPTION AND BOUNDARY EFFECTS 



Sound energy is "lost" by absorption in the water and transmission 
into the bottom. In both cases a primary variable is frequency but the 
angle of travel of the mode and the type of bottom must be taken into 
account. This paper limits discussion to a non-scattering surface, water 
volume and bottom. It also neglects energy that is transmitted into the 
bottom and then returned to the water after refraction within the 
bottom. 

Sound absorption in the water at low frequencies (under 1 KHz) is 
due primarily to relaxation of the boric ion. Standard texts [Refs. 5,6] 
define absorption in this frequency range as: 



O.IF^ 
1 +F^ 



2v40 



where: 

a = absorptive loss in dB/m, and 
F = frequency in kilohertz. 

This loss can be converted to the fractional loss, a, for a given dis- 
tance in nepers/meter: 



Loss = 20 log-jo 
ar = 20 log-]o 

= 20ar logio e 
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= 8.7 or 



a = 8.7a dB/meter, or 



= g-y nepers/meter 



2.42 



This loss rate can be introduced into the equations derived earlier 
by adding absorption as an imaginary component of the wave number; 



kg = k + ia 



2.42 



The exponential form provides the absorptive loss factor: 

^ikcR i(k+ia)R 
e ^ = e 



ikR 

= e 



“aR 



e 



2.43 



However, all of the expressions for pressure are given in terms of the 

£ p 

horizontal range, r, not actual "path" length, R. Since the loss 

factor can be written as: 



loss = e~“ 5 2.44 

Bottom interactions cause losses and phase changes which must be 
taken into account. In their discussion of propagation loss using normal 
modes ands an impedance boundary condition, Koch and Lindberg [Ref. 7] 
begin their discussion by defining the relationship between pressure and 
its derivative with respect to depth: 
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dz 



2.45 



il<2 



1-R 
1 +R 



P, 



where 

P = fluid pressure in the water at the boundary, 
dP 

— = derivative of pressure with respect to depth, 

1<2 = vertical wave number, and 
R = complex relfection (Rayleigh) coefficient. 

The Rayleigh reflection coefficient for the boundary between two 
homogeneous fluids is given [Ref. 5, Equation 6.30] by; 

Pb^b _ 

a . ^ 3.16 

Pb*^b ^ sJ-^6b 

Pw'^w sin0y 

where 

Cb = sound speed in the bottom at the interface, 

c„ = sound speed in the water at the interface, 

= density of water, 

Pb = density of the bottom, 

0^ = angle, measured from the horizontal at which the wave strikes 
the 

bottom, and 

0b = angle, measured from the horizontal, at which the wave is 
transmitted into the bottom. 

Substituting this value of R into Equation 2.45 and simplifying: 
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dz 



P. 



- il<2 



P w*^w 
Pb“^b 



sine^ 

sinGvi 



2.H7 



But 



sinGjj = '/l-cos^0i3 



2.148 



Snell's Law relates the incident and transmitted angles by: 



°b 

cos 0b = ~ 0,^ 



So, 



sin 0(3 = 





cos^0„ . 



2.49 



2.50 



The critical angle, 0 q, is defined as that grazing angle for which 
the transmitted angle is zero. This means that: 



9w _ cos(O) 
°w ^b 



and 



£b ^ 1 

c„ cos 0(3 ’ 



2.51 



where 0(, is the critical grazing angle, measured from the horizontal. 
Therefore, from Equation 2.50, 



sin0jj = 




cos^Gvj 

COS^0Q 



2.52 
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By defi-nition, 



cos 0.^ = ^ 

and 



sin = — . 



Substituting Equations 2.52 and 2.54 into Equation 2.47: 



2.53 



2.54 



dP _ . P w'^w / _ cos^6^ 

dz ~ ^ Pb^b V cos^Oq 



2.55 



In terms of wave numbers, Equation 2.53 can be used to give: 



or 



IE 

dz 



dz 



Pw^ 

Pb'^b 



i 



P w^w 
Pb*^b 




2.56 



2.57 



Equation 2.56 shows that, for any grazing angle less than critical, 
the pressure derivative will be real whereas, for grazing angles greater 
than critical, the derivative will be purely imaginary. If the pressure 
were complex, the derivative would be complex and complicated, too 
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complicated to solve in all cases with the basic program design. Therefore, 
only the real component of pressure was used in the determination of the 
boundary conditions. 

A practical discussion of the physical significance of the reflection 
coefficient and the relationship amongst it, the grazing angle and criti- 
cal angle is included in Chapter IV. 
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III. PROGRAMS 



Implementation of theory into the computer programs required two 
distinctly different methods. The program EXACT made accurate estimates 
of the normal mode horizontal wave numbers in order to compute the 
pressure as the sum of the pressure components from each of the modes. 
Absolute values of pressure are of no concern in either EXACT or EFP 
because the purpose is to compute propagation loss. It was assumed that 
forcing function, F, had a value of unity. The calculated pressure can 
then be compared to the pressure at one meter, and from that ratio, the 
propagation loss is calculable using decibel ratios. The program FFP did 
not compute horizontal wave numbers of modes; it assumed a large number 
of wave numbers equally spaced, calculated the pressure for each, and 
summed the individual contributions to give a final pressure. 

Both programs were written using complex numbers and double preci- 
sion (6M bit words), except for subroutine ’eigens’ which was restricted 
to single precision because it relied upon IMSL routine EQRTIS which was 
available only in single precision. 

Subroutines ’modes!’ and ’modes2’ provide many common functions for 
both programs EXACT and FFP. They calculate the transmitter and 
receiver modal pressure functions by stepping upward or downward from 
one of the boundaries, computing the pressure and its derivative with 
respect to depth by means of a fourth order Runge-Kutta technique 
[Ref. 7].' 
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A. NORMAL MODE PROGRAM 



This program is attached as Annex B. EXACT computes the approximate 
modal horizontal wave numbers by means of the subroutine 'eigens'. It 
then refines the horizontal wave number estimate to the desired accuracy 
using subroutine 'modesi' or 'modes2' and computes the pressure factors 
for each mode. Its aim is to compute the total pressure and, hence, the 
propagation loss at each range of interest. 

Subroutine 'eigens' uses a matrix method [Ref. 8] to approximate 
all the real modal horizontal wave numbers as well as the first evanes- 
cent mode. At the ranges of concern, it is unlikely that more than the 
first evanescent mode would be significant. The matrix of wave number 
estimates is found in a relatively short time and, although not guar- 
anteed to be accurate, it is complete; no modes are skipped. 

The next step involves finding the exact values of the horizontal 
wave number for each mode. The procedure utilizes the fact that, for 
each mode, both boundary conditions for a depth function (u or v from 
Chapter II) will be met. By setting the depth function or its derivative 
with respect to depth to a known boundary condition at the top or the 
bottom, it is possible to find a vertical wave number that will result 
in the boundary conditions being met at the the other boundary. The 
boundary with higher sound velocity will experience an exponential decay 
in the depth function. By starting at this boundary, it is a relatively 
simple matter to solve for vertical wave numbers that result in 
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sinusoidal variations that lead to solvable boundary conditions as the 
other boundary is approached. This avoids the possible problem that 
would exist if an attempt had been made to start at the boundary which 
had the lower sound velocity. In that case it is possible that instead 
of solving for an exponentially decaying function, the function may be 
exponentially growing, and no solution is possible. Figure 3.1 illus- 
trates how a poor choice of starting boundary could lead to an unsolvable 
situation. For each of the modes in the *main^ routine calls either 
^modesl* (when sound speed is greater at the top than at the bottom), or 

*modes2* otherwise. It calls with starting estimates of the modal 

horizontal wave number; the value from *eigens* the first time and an 

improved value each subsequent time. 

^Modes1\ using the best estimate for computes the modal depth 

function at each increment of depth, starting at the surface where the 

depth function is zero and its derivative with respect to depth is set to 

1.0. The integral (Equation 2.31) is calculated at each interval and its 
final value at the bottom is used in the correction term. 

At the bottom, the depth function and derivative are compared to the 
boundary conditions expected for a mode. For a rigid bottom, a deriva- 
tive of zero is expected. For an impedance bottom, the expected deriva- 
tive will be defined by Equation 2.56. The error in the derivative, tog- 
ether with the depth function and the integral is used to calculate a 
correction for the modal horizontal wave number. 

If the surface sound speed is greater at the bottom than at the top, 
subroutine ’modes2* is called. It performs the same function as *modes1* 
but starts at the bottom, where the depth function and derivative are 
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Possible Solutions to Depth Function Starting at: 
a: Boundary with Exponential (Convergent) and 
b: Boundary with Sinusoidal Changes (Possibly Divergent). 
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set to meet the boundary conditions for an impedance bottom (Equation 
2 . 56 ). A correction is made after the subroutine has stepped to the 
surface and comparison has been made with modal boundary condition 

expectations; in this case depth pressure function equal to zero. 

’Modesi' or 'modes2' will be repeatedly called until been 

satisfactorily estimated. The 'main' routine then computes the 'incom- 
plete' modal pressure from the modal depth functions, the integral, and 
the horizontal wave number, (Equation 2.31) for each mode, storing 

them for later use. It is termed 'incomplete' because there is no range 
dependency at this point. 

After the last mode has been processed by 'modes!' or ’modes2', 
range is introduced and partial pressures for each mode are calculated, 
keeping an updated sum (complex) of all the modal pressure factors. 

These pressure factors are then converted to propagation loss. It is a 
simple matter, once all the eigenvalues have been found, to compute the 
propagation loss for a series of ranges, as can be seen in Appendix 

B, Program EXACT. 

B. FAST FIELD PROGRAM 

This program was written with the intent of completing a Fast 
Fourier Program, FFP, which would utilize a Fast Fourier Transform, FFT, 
subroutine. However, the program was only completed to the point where 
pressure factors were computed for individual ranges (as in EXACT), 
rather for a large number of ranges (as is the Intent of an FFP), The 
computer program FFP, although incomplete, is Included as Appendix C. 
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It was mentioned previously that EXACT and FFP are provided many 
common functions by 'modesV and ’modes2\ There are also some basic 



differences 


in the way the 


subroutines are used 


by the two 


programs. 


In 


FFT, since 


no exact wave 


numbers 


ar e to be 


calculated. 


there is 


no 


requirement 


for subroutine 


^ eigens*. 


For the 


same reason 


, there is 


no 


requirement 


to find corrections to 


the wave 


numbers, so 


' modes 1 ' 


and 



*modes2* are used differently. In fact, both subroutines are required by 
FFP. 

Using the incremental wave number, nA5, subroutine ^ modes starts 
at the surface with the same boundary conditions as in EXACT and steps 
down, computing the pressure function and the derivative, through the 
depth of the upper of the transmitter or the receiver, and stops at the 
lower of the two. Values are required at the lower level for later use 
with the values from ’modes2’ to calculate the Wronskian, Then ’modes2^ 
starts at the bottom with the same boundary conditions as EXACT and 
steps upward until the pressure function and derivative have been calcu- 
lated for the lower level. 

The ^main’ routine then uses the pressure functions and derivatives 
to calculate a partial pressure (Equation 2.33) i'or that wave number in- 
crement. Each partial pressure is summed (integrated) so that, when 
the last wave number has been processed, the result is the total pres- 
sur e. 

As in EXACT, FFP computes the propagation loss for each range. 
Implementation of the FFT would make it unnecessary to step through 
ranges as must be done with EXACT; propagation loss would be calculated 
for as many ranges as there are FFT points. 
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IV. TESTS RUN ON THE PROGRAM 



Because it was difficult to gauge the correctness of the program or, 
more to the point, the veracity of this author^s ideas, it was considered 
prudent to check the results at each stage of development. Although the 
process was awkward and time-consuming it proved necessary and resulted 
in some unexpected benefits; the results of some of the tests provided 
graphical insight into the phenomenon of sound propagation. By starting 
with isospeed conditions, verifications were simplified. Solutions for 
horizontal wave numbers which resulted from the computer programs were 
checked against ^correct’ solutions from simple formulae for the non- 
absorption, rigid bottom system, as well as for the system that included 
absorption in the water. Once the impedance bottom was introduced, an 
analytic approach was required to decide if corrections were applied 
properly. A simple expedient to check on the program at any point in- 
volved plotting propagation loss against horizontal range for pairs of 
modes and observing the interference distance between nulls. This dis- 
tance, when compared to the theoretical distance would highlight an 
error if there were an anomaly. Another check involved the observation 
that, for a shallow channel such as that used in the model, losses are 
generally spherical at short ranges and cylindrical for far ranges. The 
smoothed propagation loss curves could be expected to lie roughly along 
curves predicted on this basis. 

The first tests involved studying the basic building blocks of the 
program. As the data were read from input files, they were printed into 
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a new file to ensure that correct information was being utilized. 
Another simple check involved printing a matrix of the wave numbers as 
they were computed. To illustrate the importance of such a basic check, 
it is worth remarking that the half-increment method required that the 
wave number be calculated for each half- increment of depth, including a 
depth of zero. Since the variation of Fortran used does not support a 
matrix with a zeroth element, a printout proved very useful in visualiz- 
ing the situation and pin-pointing a subtle error. In many cases it was 
only by printing out all of the variables after each step that errors 
could be identified and remedied. 

A. RIGID BOTTOM WITH NO ABSORPTION 

Subroutine 'eigens' was first checked for correctness by comparing 
its constant speed solution, for each mode, to a solution which was 
known to be accurate. This accurate solution was based on the concept 
that, for a constant sound speed, with rigid bottom and no absorption, 
the pressure function would be sinusoidal, satisfying both boundary con- 
ditions. This means that the horizontal wave number for each mode (n 
always odd) can be defined in terms of the wave number, k, and the water 
depth, H: 



- (^)^ 4.1 

Because the subroutine could only be run in single precision, high reso- 
lution was not expected. Table 4,1 shows that resolution increased lin- 
early with the number of depth increments. For example, for mode 2, the 
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TABLE 4.1. SUBROUTINE 'EIGENS' - COMPUTED SQUARES OF WAVE NUMBERS AND 

ERRORS RESULTING FROM VARIATION OF NU^^BER OF DEPTH INCREMENTS. 



MODE I.'UMBER 





1 


2 


3 




NUM OF VALUE 

INCRMT ERROR 


VALUE 

ERROR 


VALUE 

ERROR 


VALUE 

ERROR 


25 


0 .2144624294 
-0.0000906215 


0.1959892982 

-0.0009006298 


0.1525131932 
-0. 0032877SS3 


0.0237859429 

-0.0614532906 


50 


0.2144174816 

-0.0000456736 


0.1955415049 

-0.0004528365 


0. 1508744489 
-0.0016490440 


-0 . 0209414381 
-0.0167259096 


100 


0.2143947361 

-0.0000229281 


0.1953157087 

-0.0002270403 


0.1500510807 

-0.0008256757 


-0.0305268121 

-0.0071405356 


200 


0.2143832949 

-0.0000114870 


0. 1952023429 
-0.0001136746 


0.1496385139 

-0.0004131089 


-0.0342955159 

-0.0033718318 


400 


0.2143775572 

-0.0000057492 


0. 1951455441 
-0.0000568758 


0 . 1494320249 
-0.0002066200 


-0.0360238122 

-0.0016435355 


800 


0.2143746840 

-0.0000028760 


0.1951171158 

-0.0000284475 


0.1493287310 

-0.0001033261 


-0.0368554579 

-0.0008118S9S 


1600 


0.2143732463 

-0.0000014384 


0. 1951028945 
-0.0000142262 


0 . 1492770720 
-0.0000516671 


-0.0372637914 

-0.0004035563 


3200 


0.2143725272 

-0.0000007193 


0. 1950957821 
-0.0000071137 


0.1492512396 

-0.0000258346 


-0.0374661565 

-0.0002011912 




T 0.2143718079 


0.1950886633 


0. 1492254049 


-0.0376673477 
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error using 25 increments was .000900, and the error was halved by in- 
creasing the number of increments to 50, This trend continued through 
the run that used 3200 increments and resulted in an error of 0,000007. 
It was decided that an error of 0.00002, which resulted from the use of 
100 increments, was acceptable. The value of each horizontal wave 
number determined in subroutine ^eigens* was intended only as a starting 
point for the more accurate subroutines ^ modes'V and ^modes2\ The 
accuracy, therefore, has significance only in that, if two modes are 
closer together than single precision accuracy limitations, there would 
be a chance that a mode will be missed completely. 

Subroutine *eigens’ was also checked using sound speed profiles that 
varied linearly, both positively and negatively, with depth. Since no 
clear-cut means of assessing the error was available, no conclusions 
could be drawn, and results are published for only the positive gradient 
(Table 4.2a), In addition, by using a sound speed profile that produced 
two strong ducts, it was possible to show that ^eigens* was capable of 
distinguishing between two horizontal wave numbers which are very close 
together (Table 4.2b). The solutions are not precise but no modes are 
missed and the necessary information is provided to the subsequent sub- 
routines where the values are calculated to the required precision. 

Next, normalized modal pressures at a large number of depth incre- 
ments were calculated using ^modes2\ These data were combined in 
Figure 4,1 and show, for each mode, the relative amplitude of the modal 
pressure function. The profiles for all three real modes as well as the 
first evanescent mode are plotted together and care must be taken to 
note that they represent the modal characteristics , not the actual 
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TABLE 4.2. SUBROUTINE 'EIGENS' - SQUARES OF MODAL HORIZONTAL 

WAVE NUMBERS FOR A. POSITIVE GRADIENT, AND B. DUAL 
DUCT FOR 50 HZ SOURCE AND 500 METERS DEPTH. 



MODE 


POSITIVE 


.DOUBLE 


NUM 


GRADIEjNT 


DUCT 



1 


0.2023466395 


0 . 2094159534 


2 


0.2017806065 


0.2092275332 


3 


0.2015439681 


0. 20S35015S0 


4 


0.2009845610 


0 . 2052S2S3 12 


5 


0.2001803173 


0.2075239702 


6 


0.1991866971 


0.2065714750 


7 


0. 1980052050 


0 . 2054225531 


8 


0.1966276551 


0 . 2040741897 


9 


0.1950384284 


0.2025220985 


10 


0.1932226702 


0.2C075155S0 


11 


0.1911713327 


0. 1987873347 


12 


0.1888803823 


0. 1965925855 


13 


0.1863476142 


0 . 1941702449 


14 


0.1835693580 


0.1915113532 


15 


0.1805375873 


0 . 1S36060595 


16 


0.. 1772381208 


0.1854427397 


17 


0.1736509107 


0 . 1320079592 


13 


0.1697523042 


0 . 1TS2860324 


19 


0.1655174308 


0.1742585547 


20 


0.1609207800 


0. 1599038333 


21 


0 . 1559344080 


0.1551959531 


22 


0.1505243249 


0. 1601038096 


23 


0.1446458558 


0 . 1545894076 


24 


0.1382386306 


0 . 1436057637 


25 


0.1312213613 


0. 1420936151 


26 


0.1234851760 


0 . 1349764974 


27 


0.1148818518 


0 . 1271523591 


28 


0.1051995705 


0 . 1184328554 


29 


0.0941103985 


0 . 1087542301 


30 


0.0810451735 


0 . 0976845549 


31 


0.0648246885 


0.0347116344 


32 


0.0419336933 


0 . 0537321079 


33 


-0.0276373313 


0 . 0469792092 


34 


0 


-0 . 0198530514 
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pressure contributions. The figure shows that, for rigid bottom, no absorp- 
tion, iso-speed conditions, the profiles meet the boundary conditions at 
the top, where pressure is zero, and at the bottom, where the derivative 
is zero. The same test procedure was repeated using subroutine 'modesi' 
and identical results were obtained. 

To illustrate the relative degree to which each mode is stimulated, 
depending on the depth of receiver and transmitter, a trial run was made 
with a transmitter at depth 13 meters and receiver at various depths. 
Figure 4.2 depicts the relative strength of each mode at each depth. 

The second mode is most strongly stimulated and mode 1 is least 
strongly stimulated, as could be anticipated from Figure 4.1. 

Routine FFP does not use exact normal mode horizontal wave numbers, 
but boundary conditions must remain consistent and sinusoidal variations 
still occur at a rate defined by the vertical wave number, k^, which, in 
turn, is a function of the horizontal wave number, 5 . The results of 
subroutine 'modesi' were compared to the values of sin(k 2 Z>) and the 
results of 'modes2' were compared to the values of cos(!< 2 (H-z<) , where H 
is the water depth and z> and z< are the depths of the receiver and 
transmitter. The comparisons showed that the subroutines worked 
properly. 

At this point, the solution using calculated horizontal wave numbers 
in EXACT was examined for accuracy. It was decided that the accuracy 
criterion would be based upon the worst case which allowed for all 
errors to be cumulative. It can be seen in Equation 2.31 that an error 
in will exhibit itself directly in the square root and in the exponen- 
tial. A small error causes only a small error in the square root but it 
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DKITH (METERS) 
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Figure 4.1. Modal Pressure Function for Three Real Modes 
and EvanscenC Mode, Individually Normalized. 
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Figure 4.2. 



Pressure Contributions Normalized 
of the Most Strongly Excited Mode 



to the Maximum Value 
for Source at 13 Meters. 
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is critical in the phase component. It was categorically decided that 
the maximum allowable phase error for each mode would be 0.1 radians, 
and that the maximum range would be 100,000 meters. This limits the 
maximum error in to 10“® for all wave numbers.. 

In assessing possible real errors, a vector plot was made for each 
of several solutions. A constant adjustment was then introduced to each 
of the horizontal wave numbers, giving them an artificial error. It can 
be seen in Figure 4.3 that for short ranges, a large error is admissible 
and, even at -50,000 meters, the limit imposed by the accuracy criterion 
is well within the allowable error. 

Vector plots were also used to illustrate how a small change of 

range can cause a dramatic change in pressure. Figure 4.4a shows the 
change of propagation loss from 830 to 840 meters. Most of the change 
is involved with phase and the amplitude change is small. This can be 
seen best in Figure 4.4b where the range change is only one meter. It 
should be noted too, that these errors have been made cumulative 

whereas it is unlikely that the errors would ever be that way. One 

would therefore expect the total error to be less than that illus- 

trated. 

A plot of loss with range for a rigid bottom. Figure 4.5, showed a 
general 20 log R increase in propagation loss close to the source and a 

r 

change of 10 log R as the distance increased. Superimposed on this is a 
strong interference pattern from the three real modes. Sudden changes 
of pressure with changing range, as were noted in Figure are seen as 

rapid fluctuations or propagation loss in this plot. A decisive check on 
the general reliablity of the program to this point involved excluding 
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Figure 4.‘ 




o 




Pressure Components for Three Real Modes Illustrating 
Cumulative Error from Errors in Horizontal Wave Number 
at 50 Meters (Upper) and 50 Km (Lower). 
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Components of Pressure Illustrating Sudd 
Changes of Pressure with Range. 
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Figure 4.4. 
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Figure 4.5. Propagacion Loss for Source at 10 Meter Depth 
and Receiver at 37 Meters in 50 Meter Ocean. 
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all but two real modes and plotting the resulting propagation loss 
curves. For each pair of modes, the interference distance can be shown 
to be related to the horizontal wave numbers by: 

a - ‘•■2 

^m <^n 



where: 

R is the horizontal distance between nulls, 

and and are the real horizontal wave numbers. 

Comparison was made between the calculated interference distance and 
the observed distance for all three interference combinations and proved 
to be correct in all cases. Figure 4.6 shows the result when the first 
and third modes were used. The cycle distance was calculated to be 96.5 
meters, which agrees with the distance from the plot. 



B. RIGID BOTTOM WITH ABSORPTION 

Absorption due to water (low frequency), although normally very 
small, was made artificially high to cheek on the mechanics of the 
program. With absorption set to an arbitrary value of 0.0005 

nepers/meter, it was noted that the imaginary component of the wave 
number increased with the increasing grazing angle associated with higher 
modes: 0.000505 for mode 1, 0.000555 for mode 2, and 0.000727 for mode 

f 

3. 

To ensure that output losses were consistent with those input, other 
propagation loss curves were drawn using a variety of rates of absorp- 
tion. Figure 4.7 shows that the difference in loss at any given range is 
8.7 times the difference of absorption rates (as expected from Equation 
2.42). For example, a change from 0.0000001 nepers/meter to 0.0001 
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Figure 4.6. Interference Pattern Between Modes 1 and 3 
for Source at 10 Meters and Receiver at 
37 Meters in 50 Meters of Water for 50 Hz. 
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Figure 4.7. Propagation Loss for Absorption Rate 0.000001 Nepers/Meter 
(Upper) and 0.0001 Nepers/Meter (Lower). 
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Figure 4.8. Propagation Loss as a Function of Depth for Sources at: 
a. 1 Meter, b. 5 Meters, c. 15 Meters, d. 25 Meters, 
e. 40 Meters, and f. 45 Meters. 
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nepers/meter causes a change of propagation loss, at 10 kilometers, of 
9.0 dB. One would have expected a change of 1.0 nepers or 8.7 dB. 

A series of runs was made to illustrate how pressure varies with 

receiver depth at a set range for a variety of source depths. From the 
graphs (Figure 4.8) it is important to note that for a source or a 
receiver at the surface, the pressure will be zero ' and the propagation 
loss infinite because none of the modes is stimulated. This does not 
change when an impedance bottom is introduced, but it would change if a 
real surface were considered. 

Using these same depth profiles, it was possible to further verify 
the program by checking for reciprocity. The propagation loss is, for 
example, the same for a combination of a source at 5 meters and receiver 

at *25 meters as for a source at 25 meters and receiver at 5 meters. 

The final set of tests for a rigid bottom involved the use of the 
FFP. At the point of testing, the Fourier transform (FFT) had not yet 

been introduced so pressure factors were being calculated for individual 
ranges. Figure 4.9 shows that, for a zero absorption loss, the only hor- 
izontal wave numbers that contribute to the total are at the exact mode 
values. Small changes of the horizontal wave number cause serious 
changes in the value close to the exact mode value. 

Introduction of absorption causes the energy to be dispersed across 
the spectrum. There are small contributions from across the wave 

number spectrum and flattening at the modes. Higher absorption causes 

greater dispersion of the ^ pressured The integral, with absorption, 
can be seen to fluctuate in Figure 4.10, staying close to zero until the 

first mode is approached. The integral increases quickly, but not 
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PRESSURE 




Figure A . 9 . 



Amplitude of Sound Pressure Factor for Individual 
Wave Numbers with No Absorption. 
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PRESSURE INTEGRAL 

0.0 6.0 10.0 16.0 20.0 25.0 30.0 36.0 




Figure 4.10. Integrated Pressure with Absorption Rate 

0.00005 Nepers/Meter at Range 10000 Meters, 
Tx/Rx at 10/37 Meters, 50 Meter Ocean, 50 Hz. 
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abruptly when a mode is reached. A sharp decline is notable when the 
180 degree phase change occurs at the mode. The final value represents 
the final pressure at the receiver. 

Various checks were made to ensure proper functioning of the FFP 
program for a given range. By changing the number of wave number incre- 
ments, it was found that the maximum amplitude for each mode changed 
(Figure 4.11) but that the integrated pressure was identical when the 
number of increments was changed from 1024 to 2048. Pressures agreed 
with one another to the seventh decimal and fifth significant digit. 
Given more time, it would have been interesting to see how few samples 
could be taken, for a single range, before significant errors were intro- 
duced. 

Figure 4.12 shows the effect of changing absorption rates. However, 
the change of pressure at a given range does not correlate with the 
change of absorption rate. This indicated that the FFP routine was in- 
correct in some way and further checks were required. 

The integrated pressures from this program were then compared to 
the values obtained from the program EXACT for different absorption 
rates and ranges. As can be seen in Table 4.3, agreement was poor. 
Because of time constraints, it was decided to abandon the FFP and con- 
centrate on solutions using the EXACT method. 

NOTE: A great deal of time had been spent on FFP trying to find a solu- 

tion for a fixed range, a prerequisite to implementation of the 'FFP' 
subroutine. At the time that work with the FFT was first suspended, it 
gave inconsistent results that did not meet any of the testing criteria. 
A small change of Increment size, from 2047 to 2048 increments, caused 
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TABLE 4.3. CALCULATED PRESSURE FOR PROGRAM FFP AND EXACT. 



RANGE 


ABSORPTION 


PRESSURE 


FACTOR 


(METERS) 


(NEPERS/METER) 


FFP 


EXACT 


1 000 


0.0000001 


0.001 72867 


0.00139872 


1 000 


0.00001 


0.001 281 78 


0.001 38953 


2000 


0.0000001 


0.00001 486 


0.00139504 


2000 


0.00001 


0.0001 0503 


0.001 381 01 
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PRESSURE FACTOR 




HORIZONTAL WAVT: NUMBER 



Figure 4.11. Amplitude of Sound Pressure for Absorption 
0.0005 Nepers/Meter Using: 1024 Wavenumber 
Samples (Dashed Curve) and 2048 Samples (Solid.) 
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HORIZONTAL WAVT: NUMBER 



Figure 4.12. Integrated Pressure for Absorption Rate 

a. 0.00001 Nepers/Meter (Lower) and 

b. 0.0000001 Nepers/Meter (Upper) Range 1000 Meters. 
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outlandish changes in the integrated pressure. The changes of pressure 
resulting from a changed absorption rate did not correspond, even rem- 
otely, with the changes expected. At that point, priorities were shifted 
and it was decided to concentrate on solving the 'EXACT' program prob- 
lems at the expense of the 'FFP'. At a later date, too late to be prac- 
tically helpful, some of the FFP problems were rectified and the above 
mentioned inconsistencies were remedied. It- was possible to relate 
propagation loss to absorption rate and it was shown that the number of 
samples has a negligible effect on the final pressure. However, the 
discrepancy between the EXACT solution and the FFP solutions still exist 
and the FFP is incomplete. 
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C. IMPEDANCE BOTTOM 



The introduction of an impedance bottom to replace the perfectly 
reflecting rigid bottom made it possible to better predict the propaga- 
tion loss experienced in reality. This was done by defining the bottom 
in terms of its density and the sound speed in the bottom at the inter- 
face. After continued failure of the program to solve for modes whose 
grazing angles were far below the critical angle and the eventual reso- 
lution of that problem, a series of tests was devised to try to optimize 
the program. These tests also provided a convenient means of illustrat- 
ing how various angles of incidence associated with the different modes 
resulted in the varying phase changes and amplitudes of reflected 
energy. 

As with the rigid bottom case, it was possible to plot the modal 
pressure as a function of depth, from the reflecting surface to the imp- 
edance bottom. Figure 4.13 shows how the pressure reaches a maximum 
well above the bottom for mode 1. The grazing angle, in this case, is 
less than the critical angle and illustrates the phase advance at the 
bottom. The test was repeated for strong positive and negative gradi- 
ents but showed little change. 

The first set of tests to check on final accuracy involved using a 
critical angle of 0.3 radians and was intended to find the optimal 
number of increments. It had been determined earlier that accuracy to 
the ninth decimal, or better, was obtainable (iso-speed) for all modes 
using 200 depth incements/mode. At 100, the accuracy dropped to the 
seventh decimal and at 50, answers were accurate to only the fourth 
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DEI Til (METERS) 

60.0 <6.0 <0.0 36.0 30.0 26.0 20.0 16.0 10.0 6.0 00 




0.00 C.25 0.50 0.75 l.OO 

RILATIVT PRESSURE FUNCTION 



Figure 4.13. Relative Pressure Function for 3 Real Modes for Critical 

Angle 0.3 Radians, Isospeed Water Conditions, 50 Hz Source. 
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decimal. Table shows horizontal wave numbers that were calculated 



using 100, 200, 400 and 600 depth increments/ mode in combination with 
correction factors of 1, 2, 5 and 10. The correction factor was divided 
into the mode wave number correction term of subroutines 'modesi' and 
'modes2'to limit the size correction and thus prevent instability . 
These correction factors were found necessary for modes whose grazing 
angles were less than the critical angle (in this case only mode 1). 
Without a correction factor, the program would fail to converge to an 
acceptable solution but would become divergent. 

A major problem in deciding the best choice of variables came about 
because there was no correct answer against which to appraise the 
values. It was reasoned that 800 increments per mode and a correction 
factor of 10 would give the most accurate answer and the problem nar- 
rowed down to finding a combination that gave an acceptable accuracy 

with a reasonable number of iterations. Having decided that 200 incre- 
ments/mode was sufficient and necessary for the rigid bottom case, it 
became a question of how large the correction factor could be. The 

choice was further narrowed to a correction factor of 2 with an error of 
9.8 X 10“® and 40 iterations required and a factor of 5 with an error of 
3 X 10“® and 89 iterations. Because the error was borderline for the 

correction factor 2, the final choice was 200 intervals/ mode with a cor- 
rection factor of 5. 

Using the inputs as decided above, an effort was made to correlate 
actual accuracy with that specified as a requirement. Results for 

trials on all modes using accuracy criteria from 10~® to 10“” follow in 
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TABLE 4.4. PROGRAM 'EXACT' - REAL AND IMAGINARY COMPONENTS OF 
FOR MODES 1 TO 4 RESULTING FROM VARIATION OF NUMBER 
OF DEPTH INCREMENTS AND 'CORRECTION' FACTORS. 



COR i:; 


C ITER 




2 


3 




1 


100 


90 


0.212512S34 

0.000450164 


0. 195507949 
-0.C03339831 


C. 150035436 
-0.010172100 


C.002S67414 

0.037779947 


2 


100 


38 


0.212512965 

0.000450172 


0. 1955C7949 
-0.C03339331 


0 . 150035436 
-0.010172100 


0.002867414 

0.037779947 


5 


100 


39 


0.212513034 

0.000450172 


0. 195507949 
-0 . OC333933 1 


C. 150035436 
-0;010172100 


0. 002557414 
0.037779947 


10 


100 


185 


0.212513060 

0.000450170 


0. 195507949 
-0.003339831 


0. 150035436 
-0.010172100 


0 .002867414 
C. 037779947 


1 


200 


74 


0. 2125123S5 
0.000450165 


0. 195507947 
-0.CC3339330 


0. 150035443 
-0.010172094 


0.002S6''435 
0 . 02 77 79560 


2 


200 


40 


0 . 212512963 
0.0C04501‘:i 


0. 195507947 
-0.C03339330 


0. 150035443 
-0.010172094 


0 . 002367435 
0 . 3 / / . 6 c 0 


5 


200 


39 


0.212513C35 

O.OOC*r5Cl"2 


C. 195507947 
-0.003339830 


0. 150035443 
-0.010172094 


0 .CG2S5-'435 
0. 03:'''T9660 


10 


200 


195 


0.212513063 

0.0C0450170 


0. 195507947 
-0.003339330 


0. 150035443 
-0.010172094 


0.002567435 

0.037779660 


1 


400 


67 


0.212512855 

0.0C0450165 


0. 195507946 
-0.003339330 


0. 150035443 
-0.010172093 


0.002867437 

0.037779542 


2 


400 


40 


0.212512959 

0.0C04501':i 


0. 195507946 
-0.003339830 


0. 150035443 
-0.010172093 


0 .002367437 
0 . 03 7779642 


5 


400 


89 


0.212513033 

0.G004501''2 


0. 195507946 
-C. 003339330 


0. 150035443 
-0.010172093 


C . 002367437 
0.037779642 


10 


400 


136 


0.212513065 

0.0CC450170 


0. 195507946 
-0.003339330 


0. 150035443 
-0.010172093 


C 002357*-3'^ 
0 . 037”79642 


1 


800 


72 


0.212512685 

0.000450166 


0. 195507946 
-0.003339830 


0. 150035443 
-0.010172092 


0.002867437 

C.C37779641 


2 


SCO 


39 


0.212512959 

0.000450171 


0. 195507946 
-0.003339S30 


0. 150035443 
-0.010172092 


C.C02367437 

0.037779541 


5 


800 


90 


0.212513039 

0.000450171 


0. 195507946 
-0.003339830 


0. 150035443 
-0.010172092 


0. 00266"437 
0.037779641 


10 


800 


187 


0.212513066 

0.000450170 


0.195507946 

-0.C03339830 


0. 150035443 
-0.010172092 


0.002357437 
0 . 037779541 


EXACT 

(RIGID 


BOTT) 


0.214371320 

0.000505340 


0.195038818 

0.000555269 


0. 149226333 
0.000725949 


0.002867437 

0.037779540 
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Table 4.5. Assuming that the answers obtained using 10“'^ had the least 
error and could then serve as a standard, errors were calculated as the 
difference between the standard and each answer. It was found, in all 
cases, that the actual errors were slightly less than those specified by 
the accuracy criteria. This was to be expected since corrections were 
made to make the error less than the designated amount. 

The final check involved the problems associated with low modes 
whose grazing angles were less than the critical angle. A series of 
calculations was made to find the number of iterations required for 
different critical angles (different sound speed ratios) using a 100 
meter deep, isovelocity ocean that supported seven real modes. It was 
also intended to uncover further problems associated with large critical 
angles. At the same time, the squared values of horizontal wave numbers 
for the seven modes were calculated for each critical angle and are 
shown in Table 4.6. For angles close to the grazing angle for the first 
mode, 0.1 radians, the program failed using both subroutines 'modesi' and 
'modes2' during these runs. This was unusual because previous data had 
been collected for almost identical conditions. 

The fraction of energy reflected and the phase change represented by 
the Rayleigh coefficient, were calculated for each mode for each criti- 
cal angle and are shown in Figure 4.14 and 4.15 (upper). The graphs show 
that the loss per bounce is zero when the grazing angle is less than the 
critical angle. Mode 1 has the lowest grazing angle and it can be seen 
that the critical angle must be very low if the low modes are to suffer 
any attenuation at all. Higher modes suffer higher losses. To show the 
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TABLE 4.5. PROGRAMS 'EX^\CT' - REAL AND IMAGINARY COMPONENTS OF 
AND NUMBER OF ITERATIONS TO MEET REQUISITE ACCURACY. 



MODE NUI'IBER 





1 


2 


3 


4 . 


ACC 


REA.E 
I MAG 

ITER.ATICMS 


REAL 
I MAG 

ITER.ATIONS 


RE.AL 
I MAG 

ICE RAT I CMS 


REAL 
I MAG 

ITERJ^TIOMS 


E-5 


0.212519578303 
0 . 000446662233 
18 


0 . 195507093057 
0.003339285283 
6 


0. 150034762633 
0.010169974314 
5 


0 . 002S674’:"1359 
0.037779636753 
3 


E-7 


0.212512933248 

0.000450155917 


0. 195507946339 
0.003339329696 
0 


0. 150035442323 
0.010172092319 
7 


0.002367437052 

0.037779642325 

4 


E-a 


0.212512891707 

0.000450175339 

38 


0.195507923776 

0.003339803692 

9 


0. 150055443413 
0.010172073562 
8 


0.0C2867437C52 

0.037779642325 

4 


E-9 


0.212512888681 

0.000450176659 

44 


0.195507928795 

0,003339804619 

11 


0. 150035441633 
0.010172073282 
9 


0.002867437109 

0..C37779642316 

5 


E-10 


0.212512SS8032 

0.000450176927 

51 


0. 19550792S955 
0.003339304648 
12 


0. 150035441537 
0.010172073444 
10 


0.0C23674371C9 

0.037779642316 

5 


E-12 


0.212512387979 

0.000450176947 

58 


0, 19550792S962 
0.003339804619 
14 


0. 150035441601 
0.010172073451 
11 


0.002367437109 

0.037779642316 

5 


E-13 


0.212512887974 

0.000450176949 

70 


0. 195507928961 
0.003339804620 
16 


0. 150035441602 
0.010172073449 
13 


0.C028674371C9 

0.037779642316 

6 
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PHASE CHANGE (RADIANS) 




CRITICAL ANGLE (RADIANS) 



Figure 4.14. Phase Change/Bounce for 3 Real Modes as a Function 
of Critical Angle. 
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TABLE 4.6. 



PROGRAM 'EXACT' - REAL 
FOR IMPEDANCE BOTTOM AS 
(0 TO .25 RADIANS) . 



AND IMAGINARY COMPONENTS OF 
A FUNCTION OF CRITICAL ANGLE ’ 



MODE NUMBER 



GRIT 

ANGLE 


1 


2 


3 


4 


5 


6 


7 


RIG 


0.214751 


0.210105 


0. 20C490 


0.165134 


0.162416 


0. 128489 


0.068307 


BOT 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 


O.COOOGO 


O.OCGOOC 


.25 


0.209522 


-0.209622 


0.200611 


0. 185266 


0.152595 


0.128830 


0.070372 




-0.000004 


-0.000004 


-0.001557 


-0.002882 


-0.004518 


-0.007198 


-0.015854 


.20 


0.210226 


0.210226 


0 . 200600 


0.185258 


0.162589 


0. 128325 


0.070365 




-0.000471 


-0.000471 


-0.001785 


-0.003022 


-0.004612 


-0.007258 


-0.015882 


. 15 


0.214297 


0.210212 


0.200591 


0.185252 


0.162584 


0. 128821 


0.070360 




-0.000004 


-0.000863 


-0.001949 


-0.003126 


-0.004682 


-0.007303 


-0.0159C2 


. 10 




U1JAVAILA3LF 

WIA.VAILABLE 










.05 


0.213628 


0.213628 


0. 198860 


0.185216 


0.162580 


0.128792 


0.07CC67 




-0.000499 


-0.000499 


-0.000504 


-0.000591 


-0.002909 


-0.005782 


-0.014147 


o 

o 


0 .214649 


0.210107 


0.200526 


0.185199 


0.162529 


0. 128739 


0 . 069993 




-0.000123 


-0.000712 


-0.001612 


-0.002578 


-0.004114 


-0.006535 


-0.014487 
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extreme modal dependence of attenuation, Figure 4.15 (lower) displays 
the loss per meter for each mode for a variety of cirtical angles. 

These curves illustrate that the increased grazing angles associated 
with the higher modes result, not only in higher loss per bounce, but 

also in an increased number of bounces. For a small critical angle, say 
0.1 radians, the loss per bounce is approximately three times as great 
for mode 2 as for mode 1, but the loss per meter is about ten times as 
great (Figure 4.15). 

It was found that a change of critical angle (change of the sound 
speed ratio) changes the loss greatly. Figure 4.16 demonstrates the 
change of pressure due to mode 1 when the critical angle is changed from 
above to below the grazing angle. Changing the impedance ratio by ten 
per cent caused a corresponding increase in the loss but did not change 
the point at which the critical angle became effective. 

To study the bottom loss effects for grazing angles greater than 
critical, different runs were made for critical angles that allowed a 

different number of modes to propagate in the low loss manner (grazing 

angle less than critical angle). Figures 4.17 and 4.18 show how much 
each mode contributes toward the total pressure, depending primarily 
upon whether its grazing angle is above or below the critical angle. 
They show that for each mode, when the grazing angle is less than the 
critical angle, the propagation loss for that mode is identical to that 
for the rigid bottom case. 

As a last demonstration of the effects of the impedance bottom on 
propagation, data were compared for the isovelocity case and the strong 
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Figure 4.15* Loss Bounce (Upper) and Loss/Meccr (Lower) 
Due to Bottom Interaction. 
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DEPTH (METERS) 

60.0 46.0 40.0 D5.0 30.0 26.0 20.0 16.0 10.0 6.0 0.0 




Figure 4.16. Pressure Function for Mode 1 With Critical Angle Above 
Grazing (Crit = 0.30) and Below Grazing (crit = 0.10). 
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PROl'AGATKJN LOSS (Dli) PROPAGATION LOSS (DB) 

120.0 110.0 100. 0 00,0 (10.0 70.0 00.0 ,2o.O 110.0 100.0 00.0 eo.o 70.0 00.0 





DISTANCE (METERS) 

Figure 4 .1?. Propagation Loss for Each Mode When Critical Angie 
is Less Than Grazing Angle for All Modes (Upper) 
and For Mode 1 Only (Lower). 
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PROPAGATION LOSS (DO) PROPAGATION LOSS (DO) 

100,0 90.Q 80.0 70.0 60.0 120.0 110.0 100. 0 00.0 80.0 70.0 60 0 





DISTANCE (METERS) 



Figure 4.18. Propagation Loss Cor Each Mode When Critical Anglo 

Less Than Graziiig Angle lor Modes 1 and 2 Only (Upper) 
and for Modes 1, 2, and 3 (Lower). . 
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positive gradient case when the critical angle was 0.1 radians. It can 
be seen in Figure 4.19 that the upward refraction casued by the gradient 
lessens loss due to bottom interactions. 
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PROPAGATION LOSS (DB) 

100.0 90.0 00.0 70.0 00.0 00.0 40.0 




Figure 4.19. Propagation Loss with Impedance Bottom for 
Isovelocity Water Conditions (Dashed Lines) 
and Strong Positive Gradient (Solid Line). 
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V. CONCLUSION 



The program EXACT, which computes propaation loss by first solving 
for individual mode wave numbers, was found to give realistic values. 
The Fast Field Program, which was meant to solve for loss using an FFT 
was not completed because of complications that proved unsolvable in the 
limited time available. 

One of the problems with the solution of the program EXACT, which 
would also be a problem with FFP, is the restriction placed on the impe- 
dance bottom. All of the sound transmitted into the bottom is consi- 
dered lost to the bottom. No allowance is made for this energy to be 
refracted upward and re-transmitted into the water. This might not be 
entirely unrealistic in some cases, because bottom absorption is often 
very high. However, this simplification limits the practicality of this 
program and is considered by the author to be its main weakness. 

Another point of consideration is the number of modes involved. The 
tests run usually Involved only three real modes, two of which are lost 
to the bottom after a short range when the critical angle is 0.1 radians 
(sound speed ratio is 1.005:1). An increase of depth or frequency will 
mean an increased number of modes. There might be M, 400, or 40,000 
modes present but a large portion of them may suffer a large bottom 
loss. This also explains why trapped modes are so significant and why 
it is often sufficient to consider only a few modes in most analyses. 
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One of the intended purposes of the thesis was to compare the com- 
puter processing time for each program to solve for conditions where a 
large number of modes were present. It was reasoned that an increase 
of frequency and/or depth would require a proportionate increase in the 
number of depth samples used by the subroutines that calculated pres- 
sure as a function of depth for each horizontal wave number. This would 
mean that if the frequency were increased from 50 to 100 Hertz and the 
depth from 50 to 500 meters, there would be a 20 fold increase in the 
processor time used to do these calculations in each of EXACT and FFP. 
However, there would be a further 20 fold increase in time for EXACT 
because it would be required to calculate for 20 times as many modes. 
The FFP would not require any similar increase and would therefore have 
a decided advantage when a large number of modes were involved. It is 
possible that for cases with a small number of modes, EXACT would be 
faster, but the information was not available for comparison. 

A simplification in EXACT involves the accuracy criterion of the hor- 
izontal wave number. It was decided that errors due to accuracy limita- 
tions would be acceptable if the wave number were resolved to less than 
10~®. However, the subroutines 'eigens', 'modesi', and 'modes2' return 
the square of the wave number. For the calculated to remain within 
the limit. 



+ I Error | < ( 5^ + 1 0"®) ^ 

SO 

I Error | < 2 x 10~® 
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For the conditions considered, the correct accuracy criteria should have 
been marginally more restrictive for all three modes considered. Had 
any very small wave numbers resulted, the allowable accuracy criteria 
should have been much more restrictive. A more complete program would 
calculate each wave number to its own allowable accuracy limit using 
formula 5.1 

While it is recognized that program results would in no way alter 
the facts of nature, the correspondence between observed phenomena and 
predicted results did add to the credibility of the theory and program 
outlined in this thesis. Many of the results, as illustrated in the 
graphs, would prove useful as teaching aids to elucidate the principles 
of sound propagation in the ocean. Without reference to the development 
of normal mode theory, it is possible to illustrate the consequences of 
changing depth, range, frequency, water depth, and sound speed profiles. 
For example, the program could be used to show -the effects of surface 
duct on long range propagation. 

Finally, it is realized that the programming methods used are far 
less efficient than many readily available programs. The intention of 
this thesis was to develop an individual Normal Mode program from a the- 
oretical basis without direct reference to other, more sophisticated pro- 
grams. A more complete study would compare program results to actual 
observed data as well as to results of other, more authoritative 
predictions . 
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APPENDIX A 



MATRIX METHOD FOR FINDING MODAL HORIZONTAL WAVE NUMBERS 

The matrix method (Ref. 8) employed in the program EXACT yields a 
complete set of eigenvalues in a single pass utilizing IMSL routine 
EQRTIS. It is convenient, fast and has the advantage that all modes are 
assuredly determined, no matter how close together they are. Accuracy 
is not attempted; that is assured with the Runge-Kutta routines that are 
used subsequently. 

The Matrix solution is a means of solving for the eigenvalues of the 
one dimensional differential equation: 




where: 

(}) is the eigenfunction. 



1<2 is the vertical wave number , 

5 is the horizontal wave number, and 
z is the depth. 

It solves the equation by finding the eigenvalues of the matrix M in 
the following Matrix equation: 



miji = -(Az)^ Ei|) 



A. 2 
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This equation results from writing the differential 
difference equation at each of M depth Az apart. 

M elements ,• each element defined by: 



equation as a finite 
a column matrix with 



= ii(nAz) n = 1 , 2 N 



A. 3 



and M is a square symmetric tridiagonal matrix which allows a computation- 
ally efficient means of finding the eigenvalues. 



D, -1 0 

-1 Dj -1 

0 -1 Dj 



0 

0 



0 

0 



where: 



A.i) 



-1 Dm -1 -1 

0 - 1 Df J 



Dn = 2 - (Az)^k^(nAz) 



n = 1, 2, ,N 



A. 5 
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APPENDIX B 



COMPUTER PROGR/uM EXACT 



rlLi:;: .-KTKAM .M 



C- MAIN PROGRAM FRADS UP TO 30 PAIRS OF SOUND VELOCITY INFORMATIO.': 

C- IT ASSIGNS CALCULATED VALUES OF V/AVENUMBER FOR EACH HJCREMENTAL 

C- DEPTH AND CALCULATES THE NUMBER OF REAL EIGENVALUES AS v/ELL AS 
C- MAX VmVENUMBER AND THE DEPTH INCREMENT 

C* ABSORPTION IS ASSUMED «T0 EE CONSTANT AT ALL DEPTHS 
C’^ NON-COMPLEX EIGENVALUES ARE CALCULATED BY A MATRIX METHOD lU 

C- SUBROUTINE EIGEMS AND THESE ARE USED AS A FIRST ESTIMATE IN CL'E 

C- OF MODES SUBROUTINES 

C- PRESSURES AND THE WRONSKIAN FROM MODES ARE USED TO CALCULATE 
C- COMPLEX INDIVIDUAL MODAL PRESSURES V/HICH ARE SUMMED FOR 
C* PARTICULAR RANGES 

C- BOTTOM ABSORPTION NOT ACCOUNTED FOR 

C* 

Q-x -K -X yr yr -r -k -k * -k * -k y<: -it ■k -h * ir ic -k -k -k if -K -i- -K yr -k * -X -x y<: ^ -K -k -K -r; if -V if -k ^ yr yt -r * -k rr -k 

COMPLEX K( 10001) , P ( 10001 ), CALC ( 50 ) , ARGUE ,VAL , CCRR , INT, ?D I FF , ?TX , 
1?RX,DENOM,PPRESS(50) ,EIG(50) ,PRESS(50) ,PRESSR, IMAG , ALFA , AIG ( 50 ) , 
IFACT, REF, COMP 

REAL D ( 30 ) , KAY ( 30 ) , RATE ( 30 ) , C { 30 ) , E IG2 { 50 ) , MAG ( 50 ) , KMAX, K1 , 

1KAE( 10001 ) , KAA( ICOOl ) , PR(5C) ,PIM(50) 

LOGICAL TO? 

C 

TOP = .FALSE. 

NUM=1 

PI = 3 . 141592654 

CMAX=0 

CMIN=200 

FREQ = so 

DEPTHR = 10 

DEPTHT =37 

RMIN =20 

INC = 400 

JUMP = INC/100 

ERR.MAX = .00000001 

CORN = 2 

CRITA = .15 

R.ATIO = 2.0 

12 RE.YD (3,’^) D(NUM) ,C(NUM) 

C 

Q-:vky^y:kkkkkkkk-/<kkkkkr‘ki-k'tr-:^kr(^Tr’Kkk-kkkk-krfir-*-7^-r-f<k-kk-kk7f:k-r-kkk-kk-kk-K':^-k-x-x-^-^ rr 

C- EACH SOUND VELOCITY/DEFTH PAIR HAS A V/AVENUMBER CALCUL.^TED 

Qi( k k k k k k k ic k k k k k k k k k k ir if k k k ^ k k ir ^ ir k k -:r k k k k k k k k k -k k k k k k -k k k ■x k k k -K -k -K k k ^ k - k y- 

c 

IF (D(NUM) .LT. 0) GOTO 14 
KAY(NUM) =2 - PI - FREQ / C(NUM) 

IF (C(NUM) .GT.CMAX) THEN 
CMAX=C ( NUM ) 

END IF 

IF ( C ( MUM ) . LT . CM I M ) THEM 
KMAX=KAY(NUM) 

CM IN = C(NUM) 

END I F 

IF (NUM .EQ. 1) GOTO 13 
C 

R.RTE (NUM-1) = (KAY(NUM) - KAY(MUM-l)) / ( D ( NUM ) -D ( MUM- 1 ) ) 
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PRINT 103 , D ( NUI-> 1 ) , C ( MUN- 1 ) , KAY ( NUN - 1 ) , RATE ( NUM- 1 ) 

13 NUM = NUM ^ 1 
GOTO 12 

14 BOTTOMED ( NUN- 1 ) 

N = IFIX( (4-^B0TT0N-FRE0/CNA:< - l)/2.0) 

VJRITE(4, 6389) M 

6339 FORMAT (//’ THERE ARE '.16,' REAL MODES'/) 

M = 1 + N 

DELZ = BOTTOM/ ( I NC*N) 

PRINT 103,D(NUM-1) .C(I'aJM-l) 

V/RITE (4, 103 )D(NUM-1 ) ,C(NUM-1 ) . KAY(NUM-1 ) 

IF (C(NUM-l) .LT.C(l) ) TOP = .TRUE. 

IF (M.GT.50) THEN 
WRITE (4, 100) N 

100 FORM.AT('TOO MANY MODES FOR THIS PROGRAM!! M(REAL) = ',14) 

GOTO 11 
END IF 

X •rxx'A-xxxxxxxxxxxxxxxxxxxxx^ xxxxx'xxxxxx 

* 

THE NEXT SECTION CALCULATES THE EXACT VAAl’ENUMBERS FOR THE SU3- 

* ROUTINES El GENS, MODES 1 AND M0DES2 

* THE REAL COMPONENTS FOR MODES 1 AND MODES 2 ARE CALCULATED FOR 

" INC*N + 1 DEPTHS STARTING WITH DEPTH 1 AT THE SURF.ACE AliD INC- 

- AT THE BOTTOM 

" AT REGULAR DEPTH INCREMENTS A REAL VALUE CF K IS TAKEN FOR 

USE IN SUBROUTINE EIGENS THE NUMBER OF SAMPLES IS ICO 



KOUNT = C 
NUM = 1 

DO 18 IG = 1. INC-M+1 
DP = (IG-1) - DELZ 

K A A ( I G ) = KAY ( NUM ) RATE'. ( I iUM ) - ( DP - D ( I :UM ) ) 

IT IS POSSIBLE TO INSERT A LOOP TO CHOOSE ONLY A FRACTIC:: OF 
THE DEPTH SAMPLES OF K( ) A TEST RUN WAS MADE TO CC.MPARE THE 
VALUES OBTAINED FOR VARIOUS 'JUMPS' 

rx'x*'A’:fe':^xxxxxxxx-^x-^x:txxxxxxxx*xxxxx‘^x:A:*X'Axr.'^xxxxx'r 

IF (INC.LE.200) JUMP = 2 
I F ( MOD (IG-1, JUMP ) . EQ , 0 ) THEN 

KOUNT = KOUNT + 1 
IF(IG.EQ.l) THEN 
KOUNT = 0 
GOTO 886 
END IF 

KAE( KOUNT) = KA.A(IG) 

886 END IF 

IF (DP.GE.D(NUM + 1)) NUM = NUM + 1 
18 CONTINUE 
C 

CALL El GEMS (BOTTOM. KOUNT, M. KAE, EIG2 , RMIM) 

C 
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■K TT 



I-:XACT 



'I 



Q -A- -r ^ X -V X X -X :v w rf X r r x 

C- ABSORPTION IN NEPERS/'MTTSR CALCULATED AS FUNCTIOrj OE EREQ 



FREK = FREQ/ 1000 • 

FREER = FREK- ^2 

ALFAW = . 00001 *FREEK/( l^FREER) 

C 

Q S: A X '.r * 7c X X Vr ★ T- T>r X X X x ■x X X ★ X X w -^r X- ★ TT x >• tV X x x -r x x x x x x 'Hr x ••■.• x x x x x x •*■ x -v x x x x x x x x x x x x 



C- NEXT FOLLOV/S A LOOP WHICH CAN EE USED TO DETERMINE THE PRESSURE 

C* OR ATTENUATION AT ALL DEPTHS FROM TO? TO BOTTOM 

QHrxxxx:txxxxHcHrxxx^HrxxHrHcx-xxxxTrHrxxxxx*x‘*-xxxxxxxxx'xHrxxxxxx*Hrxxx 

c 

DO 54 IDR = 37,37 
DEPTHT = 1.0-^ IDR 
DO 55 I =1,M 
CRIT=CRITA 
NC = 1 
ALFAB =0.0 
ALFAT = ALFAW ALFAB 
ALFA = CMPLX(0.0, ALFAT) 

5500 DO 598 IK = i,INC-N+l 

K(IK) = KAA(IK) + ALFA 
598 COMTIN^UE 

ARGUE=K ( 20)--2-(PI-(I-.5) /BOTTOM ) * ^ 2 
CALC(I) = CSQRT( ARGUE) 

0**xx**HrHc*^xxxx:^r;TxHr>/Hrxx*HrxxH:x*'X'x^*xxxHrxxxxrVxxxxWxHrx + x* 

C^ CALC IS USED TO EVALUATE THE CORRECTNESS OF THE PROGRAM 

C^ BY INPUTTING AN ISOVELOCITY PROFILE IT IS POSSIBLE TO CALCULATE 

C* THE CORRECT HORIZONTAL L\AVENUM3ER FOR EACH MODE AND THEN COMPARE 

C- THE RESULTS OF THE PROGRAM TO THIS 'CORRECT' ANSWER 

Qxxx**7--:V'xxxxxxxr‘:**XX*xxx'xx-xxx7.xx:<rxx**xxxx-^x xxxxx 

VAL = EIG2 ( I ) 

CORRIJ = 1 ■ 



IF (TCP) THEN 

Qxxitxxxxx'*‘xxx*xxxxxxxrVxxxxxxxxxxxxxxx 



C- IF THE SOUND VELOCITY IS GREATER AT THE TO? THAN AT THE EOITCM 
THEN SUBROUTINE MCDESl IS USED; OTHERWISE MCDES2 IS USED 
AFTER EACH RUN OF EITHER SUBROUTINE, THE CORRECTIOIi MADE TO THE 
C- SQUARE OF THE HORIZONTAL V.’AVELENGTH IS COMPARED TO A ERROR 
C* CRITERION; IF IT EXCEEDS THE SET LIMIT THE SUEROUTII^E IS CALLED 

C’^ AGAIN. IF THE CORRECTION IS LESS THAN THE LIMIT, THE PROGR.A.M 

C- PROGRESSES AND SOLVES FOR THE PRESSURE FOR TK.AT MODE 

Qif * rr ■k * i: if iv 7C y: i: -k * -k <■ -k -k k -k ir -ir -k k ir -k -k ir k; -k -k i: T- 7^- k if -k- re 'fr -k -k if ir kr :r -k -k -k -k k ^ ^ -r k: ir -k ■*- -k -fc -k -f: 

c- 

180 CALL MODES 1 (VAL, DELZ, M, K, P,CORR, I NT, INC, RATIO, CRIT, CORRN, I 

NC = NC + 1 

C IF(CABS(VAL) .GT.CABS(K( INC-N^ 1 ) .^’^ 2 ) ) THEN 

C CRIT =.5-CRIT 

C CORRN = CORN 

C IF (CRIT. LT. .005 ) GO TO 66 

C VAL = EIG2 ( I ) 

C GOTO 180 

C END IF 

IF(CABS(CORR) .GT.ERRHAX) GOTO 180 
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^XA'JT 



T : LE : 



I:vRT:vAN 



A1 



IF (CRIT. LT.CRITA) THEN 
CRIT = CRITA 
GOTO 180 
END I F 
GOTO 57 
WRITE(4, 570) I 

I F { VAL . EQ . 0 ) THEN 
ELSE 

190 CALL MODES2(VAL,DEL2,M,K,?,CORR, INT, INC , RATIO, CRIT, CCERN, I ) 

NC = NC + 1 

I F( CABS (VAL) .GT.CABS(H( ) - -2 ) )THEN 

CRIT =.5-CRIT 
CORRN = CORN 

IF (CRIT. LT. .005) GO TO 56 
VAL = EIG2(I) 

GOTO 190 
ENDIF 

IF( CAB S ( CORR ) .GT. ERRN.AX) GOTC 190 
IF (CRIT. LT.CRITA) THEN 
CRIT = CRITA 
GOTO 190 
ENDIF 
GOTO 57 

56 WRITE(4, 570) I 

570 FORMAT ( ' NO CONVERGENCE ! ! ! ! I FI 13 THE ERROR' .16) 

5 7 ENT> IF 

EIG( I ) = CSQRT( VAL) 

IDTl = IFIX(DEPTHT/(DEL2-2 ) )-2 

IDT2 = IDTl + 2 

PDIFF = P(IDT2) ~ P(lDTl) 

DDIFF = DEPTHT - (IDTl) -DELE 

PTX = P(IDTl) PDIFF-DDIFF/(DELZ-2 ; 

IDRl = IFIX(DEPTKR/(DELZ-2) ).-2 

IDR2 = IDRl + 2 

PDIFF = P(IDR2) - P(IDRl) 

DDIFF = DEPT.HR - (IDRl) -DELE 

PRX = ?(IDR1) + PDIFF-DDIEF. (DELE-:) 

AIG(I) = CMPLX(KEAL(EIG( I ) ) , A53(.-.:XAG(EIG( I ) ) ) ) 

DENOM = INT*CSQRT(AIG( I )-?I/2 ) 

PPRESS( I )=PTX-PR.X/DENOM 

C FOLLOWING INSTRUCTIONS ONLY FOR CALC OF RAYLEIGH COEFFICIENT 
C DONE FOR RANGE OF CRITICAL ANGLES .AND FOR DIFFERENT IXPEDAI.'C 

55 CONTINUE 
C 
C 

C- THE FOLLOWING LOOP CALCULATES THE ATTENU.ATION OF THE SIGNA 

C- DIFFERENT RA.NGES BY SUMMING THE COMPLEX PRESSURE DUE TO TH 

C- INDIVIDUAL MODES 

C* CAM BE SET UP TO PRODUCE ATTENUATION AT REGULAR R.ANGE INTERVALS 

C^ OR AT IRREGULAR INTERVALS SUCH .AS EVERY 10 METERS OUT TO ICO 

METERS, THEN EVERY 100 METERS TO 1000 A.ND SO ON. 
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IMAG=CM?LX(0, 1 ) 

DO 560 IRKG =1,1 

DO 561 NRMG = 300,3000 
R=( 10--IRNG)-MRXG 
C R = 50000 

FRESSR = 0 
STRONG = 0 
DO 562 I = 1,N 

PRDSS( I )=PPRESS( I )*CEXP( IMAG-EIG( I ) - R ) /SQRT ( R ) 
MAG(I)= CABS(PRESS( I ) ) 

FRESSR= PRESSR + ?RESS(I) 

PR( I )=REAL(PRESSR) 

PIM( I ) = AIMAG(FRESSR) 

EASE = ATAiN(AIMAG(PRESSR)/REAL(PRES3R) ) 

STRONG=CABS ( PRESSR ) 

C WRITE ('4,3434)1, PRESS ( I ) , PFRESS ( I ) 

3434 FORMAT ( I6,4E15.7) 

C PRINT 3487, MAG( I ), STRONG 

3487 FORMAT(2F15.5) 

FACT = CSQRT(K( lNC-N+1 ) ’^-2 - ( E IG ( I ) /COS ( CR ITA 2 ) 

REF = (RATIO - FACT)/(RATIO FACT) 

FRACT = CABS (REF) 

COMP = CSQRT(K( INC-N+ 1 ) ^’ -2 -EIG ( I ) - -2 )/(2-K( I NC^N+ 1 ) - BOTTOM 
QLOS = (1 -FRACT) -CABS (COM?) 

PHAS = ATAM( AIMAG(RSF)/REAL(REF) ) 

THETA = ARCOS(EIG( I )/K( INC-N+1 ) ) 

C 

562 CONTINUE 

510 CONTINUE 

/\TTEN = - 20- AL0G10( STRONG) 

561 CONTINUE 

560 CONTINUE 

54 CONTINUE 
1 1 STOP 
END 



SUBROUTINE MODESl ( W\L, DELE , N , K, ? , CORR, IN’T, INC, RATIO, CR IT, CORRN 
COMPLEX K( 10001 ) , P( 10001 ) , VAL,CORR, INT, FI , Al, 31, F2 , A2 , B2 , F3 , A3 
1B3, F4, B, A, IMAG, AP,DENO 
C 

C- THIS SUBROUTINE USES A FOURTH ORDER RUMGE KUTTA TECHNIQUE 
C* TO CALCULATE PRESSURE AMD ITS DERIVATIVE STARTING AT THE SURF 

C- AND STEPPING DOWN A DISTANCE 2 TIMES THE DEPTH INCREMENT UNTIL 

C- THE BOTTOM IS REACHED 

C- CALC AT EACH L IS IN FACT FOR THE L+1 STEP 

C- AT THE BOTTOM THE CALCULATED PRESSURE, DERIVATIVE AND INTEGRAL 
C- ARE USED TO DETERMINE THE CORRECTION TO BE APPLIED TO THE 

C- SQUARE OF THE BEST GUESSS OF THE NORMAL MODE H0RI2CNTAL WAVE- 



, I) 



ACE 
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::'jr*:3ER the correction is applied and the result 

TO THE NAIN PROGRAM 

IN ADDITION THE PRESSURE AT ALL DEPTH INTERVALS IS 
THIS IS USED IN THE CALCULATION OE ATTENUATION 
--SHOULD SAVE ONLY THOSE PRESSURES REQUIRED NOT; ALL 



PI = 3 . 141592654 
IMT =0.0 
E = 0.0 



C 



c 



c 



c 



P(l) = B 

DO 18 L = 1, I.NC*N - 1,2 



D? = DELZ*(L+1) 

FI = -B - (K(L)"^*2 - VAL) 

A1 = A + DELZ’^Fl 
B1 = B + DELZ-A 

F2 = -B1*(K(L+1)'‘*2 - VAL) 

A2 = A + DELZ - F2 
B2 = B DELZ -Al 

F3 = -B2*(R(L+l)--2 - VAL) 

A3 = A + 2-DELZ-F3 
B3 = B ^ 2-DELZ-A2 

F4 = -B3*(K(L+2)*^2 - VAL) 

B = B + 2’=^DELZ'^(A + Z’^Al + 2-A2 + A3)/6 
P(L+1) = B 

A = A + 2-DELZ^(Fl 2-F2 + 2’^F3 + F4)/6 
INT = INT + B’^-2^'DELZ-2 



IS PASS 
PASSED ; 



'•N^l 



c 

18 CONTINUE 

I MAG = CMPLX(0.G, 1.0) 

FILL = REAL(K( INC-N-^1) )*-2 
DENG = CSQRT((FILL - VAL)/FILL) 

AP = CSQRT ( F I LL-VAL/COS ( GRIT ) ^-^ 2 ) -E- I MAG/RAT I O-DENO 
CORR = B-(A+AP)/INT 
VAL = VAL - CORR/CORRN 
21 RETURN 
END 



SUBROUTINE MODES2 ( VAL , DELZ , N , K, ?,CORR, INT, INC, R.ATIC , GRIT, C 
COMPLEX K( 10001) , P( 10001 ) , VAL, CORR, INT, FI , A1 , B1 , F2 , A2 , B2 , 
1F3, A3,B3, F4,B, A, IM.AG,DENO 
C 

C- THIS SUBROUTINE USES A FOURTH ORDER RUNGE KUTTA TECHNIQUE 
C* TO CALCULATE PRESSURE AND ITS DERIVATIVE STARTING AT THE B 

C* AND STEPPING UP A DISTANCE OF 2 TIMES THE DEPTH INCREMENT U 

C* THE SURFACE IS REACHED 

C* CALC AT EACH L IS IN FACT FOR THE L- 1 STEP 

C* AT THE TOP, THE CALCULATED PRESSURE, DERIV.ATIVE AND INTEGR.A 

C- ARE USED TO DETERMINE THE CORRECTION TO BE APPLIED TO THE 
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C- S>^UARE OF THE BEST GUESSS OF THE NORMAL MODE HORIZONTAL 

C- NUMBER THE CORRECTION IS APPLIED AND THE RESULT IS PAS 

Cx TO THE MAIN PROGRAM 

C- IN* ADDITION THE PRESSURE AT ALL DEPTH INTERVAl^S IS PASSED 

C- THIS IS USED IN THE CALCULATION OF ATTENUATION 



PI = 3 . 141592654 
I NT = 0 

LFLOOR = INC-N + 1 
B = 1. 

P(INC’^N) = B 

IMAG = CMPLX(0.0, 1.0) 

FILL = REAL(K( LFLOOR) 

DENO = CSQRT((FILL - VAL)/FILL) 

A = CSQRT(FILL-VAL/C03(CRIT)--'^2)-B-IMAG/RATI0/DEN0 
A = 0. 



IF 


(I 


. Cl 


0 

II 

< 




DO 


28 


L 


= 1,INC-N-1,2 






DP 


= 


DEL2-^(LFLOCR-L-2) 






FI 


= 


-B * (K(LFLOOR-L+l)* 


-2 - VAL) 




A1 


= 


A + DELZ-Fl 






B1 


= 


B + DELZ-A 






F2 


= 


-Bl- (R(LFL00a-L)*-2 


- VAL) 




A2 


= 


A + DELZ-F2 






B2 


= 


B + DEL2-A1 






F3 




-B2-(R(LFLOOR-L)*-2 


- VAL) 




A3 


= 


A + 2--DELZ-F3 






33 


- 


B - S^DELZ-A2 





F4 = -53- (K(LFLOOR-L-l)^-2 - *VAL ) 

B = B + 2*DELS-(A + 2*A1 + 2- A2 + AS )/6 
?(LFLCOR-L-2) = B 

A = A + 2-DELZ-(Fl ^ 2-F2 + 2-F3 + F4)/^6 
INT = INT + B---2-DELZ--2 

23 CONTINUE 

CORR = 3--A/INT 
VAL = 7AL + CORR/CORRN 
21 RETURN 
END 



SUBROUT I NE E I GENS ( BOTTOM , KOUNT , M , KAE , E I G2 , RM I N ) 

DIMENSION E( 10001) ,D( 10001 ) ,AH2( 10001) ,EIG2( 50) , EIG(50) 
REAL KAE( 10001) 

C PRINT 4356, BOTTOM, ROUNT, M, RMIN 

it it if if if ir -k 'k -h -ir i( -k if ji -f; ic if -K i: ^ -fr -k ^ "h -k 7; i: i< if -k 'K -k -K 'r’ if it ie k -k -k -ir -K -K -k i: -k -k 

c- 

C* COMPLEX NOT USED IN THIS SUBROUTINE BECAUSE IT CALLS IMSL 
C- SUB EQRTIS AMD IT IS IMPRACTICAL TO MODIFY THAT FOR CCMPL: 

C- THE REAL EIGENVALUES ARE THEREFORE CALCULATED AND USED IN 

C- MODESl AND MODES2 V/HERE COMPLEX EIGENVALUES ARE USED 
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JSED IIs THE EFT PROGRAM 



DZ = BOTTOM /(ROUHT-i) 

DZ2 = DZ’^-DZ 

DO 100 I=l,KOUMT 
AK2(I) = KAE(I)--2 . 

D( I ) = -AK2 ( I ) -DZ2 + 2.0 
E( I ) = 1.0 
100 CONTINUE 

D ( KCUNT ) = D ( KOuMT ) - 1 . 0 
E(l) = 0.0 
ISW = 0 
lER = 0 

CALL EQRTIS ( D , E, KOUNTM’L 1 3V;, lER ) 



DO 6000 J = 1 . M 
EIG2(J) = -D(J)/DZ2 



ALTHOUGH THE FORMULA IS EXACT THE SUBROUTINE RUNS IN SINGLE 
PRECISION SO ONE MUST BE CAREFUL 

K -x yc ■x TV -K -k ^ -r ir -k -k -X f i; 'IT Tf 7,':\'X'*'vx‘^xxxr;'Txx4iVx'xr;txrr*'x:rwxx7rTC':rT«rx:^>r'X'X>'xrxxx'X‘:rrr'”-r-r-*' 

EIG2EX = AK2(20) - ( 3 . 1 415926535* ( 2- J- 1 ) /( 2 . O^BOTTOM ) ) ' - 2 
6000 CONTINUE 
C 

9000 FOR.M.AT (1 5, 3 FI 6. 10) 

8000 RETURN 
END 
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APPENDIX C 



COMPUTER PROGRAM FFP 



FILE: FFP 



FOPT 



COMPLEX WRON, U, DU , V, DV , VI , DVl , K( 25000 ) , FN , E I OVAL , VAL , I MAG , FVMC , 
lAREA, SIMPLE , SOLN, VSE , VEW , VRTV/N , WRCMR 
REAL KMAX,D( 30) , KAY (30) , RATE ( 30) , C(3G) 

C 

MUM=1 

PI = 3 . 14159265 
CMAX=0 
CMIN=2000 
C 

FREQ = 50 
DEPTKT =10 
DEPTHR = 37 
RMIN = 40 

C INC , THE NUMBER OF INCREMENTS PER REAL MODE, MUST BE EVEN) 

INC = 200 

12 READ (3,^) D(HUM) ,C(NUM) 

C 

IF (D(NUM) .LT. 0) GOTO 14 
KAY(NUM) = 2 * PI * FREQ / C(NUM) 

IF ( C ( NUM ) . GT . CMAX ) THEN 
CMAX=C(NUM) 

DMAX=D(NUM) 

END IF 

IF ( C ( NUM ) . LT . CM I M ) THEN 
KMAX=KAY ( NUM ) 

DKMAX=D ( NUM ) 

CM IN = C(NUM) 

END IF 

IF (NUM .EQ. 1) GOTO 13 
C 

RATE (NUM-1) = (KAY (MUM) - KAY (MUM-1)) / ( D ( NUM ) -D ( MUM- 1 ) ) 

C 

13 NUM = NUM + 1 
GOTO 12 

14 BOTTCM=D(NUM-l ) 

N = IFIX(2-BOTTOM-^FREQ/CMAX+.5) ■ 

INCTOT =■ INC’^N 

IF ( INCTOT. GT. 18000) INCTOT = 18000 
DELZ = BOTTOM/ INCTOT 
IF (N.GT.50) THEN 
WRITE (4, 100) N 

100 FORMAT(’TOO MANY MOLES FOR THIS PROGRAM!! N(REAL) = A 14) 

GOTO 11 
END IF 



VALUE OF ABSORPTION DEPENDEIJT ON FREQUENCY ONLY-NO TEMP DEPENDENCE) 

FFREQ = FREQ/1 COO 

ALFAl = .0000092/(.7 + FFREQ**2) 

ALFA2 = .0046/(6000 + FFREQ- ’^2) 



OOOOOOO O o OQ 



c j r 



A1 



F..RTR.-.N 



ALFA3 = . 000000046 

ALFA = (ALFAl ^ ALFA2 ^ ALFA3 ) * FFR£Q-’^2 
ALFA = .00001 
ALFA = 0 

EXAGGERATE ALFA TO SEE IF NUMBER OF ITERATIONS IMPCRTAIiT 
MUM = 1 

DO 18 IG = 1,INC-M + 1 
DP = ( IG-1) DELZ 

K ( I G ) = KA V ( ! lUM ) + RATE ( NUM ) - ( DP - D ( NUN ) ) 

IF (DP.GE.D(NUM +1)) NUM = NUM + 1 
18 CONTINUE 

IF (DEPTHT.GT.DEPTHR) THEN 
DUP = DEPTHR 
DLV/ = DEPTHT 
ELSE 

DLW = DEPTHR 
DUP = DEPTHT 
END IF 

lUl = IFIX(DUP/(DELZ-2) )-2 
IU2 = lUl + 2 

I LI = IFIX(DLV;/(DEL2^2) ) -2 
IL2 = ILl + 2 



ALFA IS CALCULATED FOR ONE FREQUENCY AND IS PART OF VAL(CCM?LEX) 

★ xx*xxx*xx*xxx***x*-'»rxxxxiVxx*xxx')<rxx + xxxxxxxxxxxxxxxxxxxT*rx*xxx 

AREA =0.0 
IFFP = 2047 
DO 55 I =1, IFF? 

C 

0XXXTV'JtlfrrV>tx^XXXxxxxxxxxxxxxxxxT<rxx + S:*';Vxxx:V*xx;rw V'x 

c 

C DURING THIS PHASE OF THE PROGRAM MODES2 PROVIDES VALUES FOR 

C PRESSURE AND DERIV:\TIVE AT BOTH THE TRANSMITTER AND TH:E RECEIV'ER 

C STARTING V/ITH THE INITIAL CONDITIONS AT THE SURFACE 

C (IE PRESSURE = 0 AND DERIV IS AN ARBITRARY VALUE) 

C MODES 2 PROVIDES VALUES FOR PRESSURE AND DERIVATIVE FOR THE 

C LOWER POINT ONLY STARTING WITH THE INITIAL CONDITIONS AT THE 

C BOTTOM (IE DERIVATIVE = 0 AND PRESSURE IS AM ARBITRARY VALUE 

C A VALUE FM REPRESENTS THE CONTRIBUTION OF EACH ' I MCREMEN'T ' 

C AND THE AREA UNDER THE FN CURVE REPRESENTS THE TOTAL ’PRESSURE 

C 

Qxxxxx':^***xxx:V-*****x + xxxx*Ttxxxxxxxx*-^x 

c 

R = 500 

WVNMCH = 1 . 1-KMAX/IFFP 

V/VNM = I* WVNMCH 

EIGVAL = CMPLX(V;VNM, ALFA) 

VAL = EIGVAL*-2 

CALL M0DES1( VAL, DELZ, INCTOT, X, U, DU, lUl, I U2 , DU? , DLW , VI , DVl , 
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FOixTRAM 



FI LZ ; FFP 



1 IL1.IL2) 

CALL MODZS2( VAL,DELZ, INCTOT, K, V, DV, I LI, IL2,DLV/, INC) 

WRON - V^DVl-Vl-DV 
878 FORMAT (10E3.1) 

7( -X i; -k -h -r -k if It 'k -k -k -k -k ic -k -k ic it -k ir -k ic ic -k ir -k :ir -.r -r if -r -k -k ■k It -k -K v ^ -K -k -T 'X yr -r -r -k -r -K -r 

C CALCULATE PARTIAL FIELD INTEGRAL FOR PARTICULAR VALUE OF 

C WAVENUMBER COMPARE TO THAT CALCED BY PROGRAM 

■K -k yr ^ -k ie •k yr -k Tf Tt -k -k K -k "X TT ir TV "K Tt -k ie -k -k i( •x TV ic TV "ir k 'K -r -r- rr yr Tc -k Tc Tf •k Tc k "k -k -k -K -k Tr -X 

C VRTV;N = CSQRT(K(20)*-2 - VAL ) 

C YEW = SIN( VRTWK-DU?) 

C VEE = COS( VRTWN- (BOTTOM-DLW) ) 

C WRONK = VRTWM'^COSCVRTV/N’" BOTTOM) 

C SOLN = YEW * VEE / WRONK 

C SIMPLE = U^V/WRON 

FN = U-V-CSORT(EIGVAL)/WRON 
I MAG = CMPLX(0,1) 

FUNC = FPF'^CEX?( IMAG^EIGVAL-R)^SQRT(2/(PI-R) ) 

CONT = CABS (FUNC) 

AREA = AREA + FUNC*V;VNMCH 
TOTAL = CABS (AREA) 

IF (MOD(I,5) .EQ.O) WRITE (4,190) EIGVAL, CONT , TOTAL 
IF(MOD(I,100) .EQ.O) PRINT 190, E I GVAL, CONT , TOTAL 
55 CONTINUE 

C 

11 STOP 

END 

Q x'l^rxwxTr-K-Sr'ic'r'r'T. r'Xir-jr - 7- 

C MODESl PROVIDES PRESS AND DERIV AT BOTH UPPER AND LOWER 

C INPUTS 

C VAL - THE INCREMENTAL HORIZONTAL W’AVEMUM - U? TO 8192 VALU 

C DELZ- DEPTH INCREMENT FN OF DEPTH, NUMB OF MODES, AND irUMB 

C OF INCREMENTS / MODE 

C N - NUMBER OF MODES 

C K - HORIZONTAL V/AVEirUMBER (COMPLEX) AT EACH DEPTH INCKEMLN'T 

C DUP - DEPTH OF UPPER OF TXER AMD KXER 

C DLW - DEPTH OF LOWER OF TXER RXER 

C IU1/IU2 INCREMENT NTjMBER FOR IMC ABOVE AND BELOW 'UPPER' 

C IL1/IL2 INCREMENT NUMBER FOR INC ABOVE AND BELOW 'LOWER' 

C INC - NUMBER OF INCREMENTS OF DEPTH 

C OUTPUTS 

C U1 - 'PRESSURE' AT UPPER 

C DUl - 'DERIVATIVE AT UPPER 

C VI - 'PRESSURE' AT LOWER 

C DVl - 'DERIVATIVE' AT LOWER 

C 

c 

SUBROUT I NE MODES 1 ( VAL , DELZ , I NCTOT , K , U , DU , I U1 , I U2 , DUP , DLW , V 1 , DV 1 , 
IILI, IL2 ) 

COMPLEX K( 25000 ), P(2 ) , DP ( 2 ) , DPI ( 2 ) , P 1 ( 2 ) , FI , A1 , B1 , F2 , A2 , B2 , 
1F3,A3,B3, F4,A,B,PDIFF,P1DIFF,QDIFF,Q1DIFF,U, V1,DU,DV1,VAL 
C 

B = 0.0 
A = .00001 
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III w 



o o o o o 



FFP 




FORTR.'.:: A’ 




DO 18 


L 


= l,I.NCTOT ^ 1, 


2 


FI 


= 


-E - (:;(L)^^2 - 


V A L ) 


A1 


= 


A + DELE'' FI 




B1 




B > DELZ -A 




F2 


= 


-B1-(:<(L-1)^''"2 


- VAL) 


A2 


= 


A + DELZ * F2 




B2 


= 


B + DELZ ’^•Al 




F3 


— 


-32’^(K(L+l)-^2 


- VAL) 


A3 


= 


A + 2-D£LZ*F3 




B3 


= 


B + 2*DELZ’^.A2 




F4 




-B3*(K(L+2)**2 


- VAL) 



B = B + 2*DELZ-(A + 2-Al + 2-A2 - A3}/6 
A = A + 2-DELZ-(Fl 2-^F2 + 2-F3 ^ F4)/6 

IF( (L+1) .FQ. lUl) THEN 
P(l) = B 
DP(1) = A 
END IF 

IF ( (L^l) .EQ. IU2) THEM 
?( 2 ) = 5 
DP(2) = A 
END I F 

IF ( (L+1) .EQ. ILl) THE-N 
Pl(l) = B 
DPl(l) = A 
END I F 

IF ( (L+1) .EO. IL2) THEN 
Pl(2 ) = B 
DPI (2) = A 
END IF 

18 CONTINUE 

FDIFF = P(2) - ?(1) 

PIDIFF = ?1(2) - Pl(l) 

QDIFF = D?(2) - DP(1) 

QiDIFF = D?l(2) - DPl(l) 

DDIFI = DUP - (lUl)-DELZ 
DlDIFt = DLW - (ILl)-DELZ 
U = P(l) + FDIFF-DDIFF/(DELZ-2) 

VI = Pl(l) + PlDIFF-DlDIFF/(DELZ"-2) 

DU = DP(1) + CDIFF*DDIFF/(DELZ-2 ) 

DVl = DPl(l) + QlDIFF-DlDIFF/(DELZ-2 ) 

21 RETURN 
END 



•K-)eyi:-tr'kyffr-ic-k-^x''-^'r-r-ir:Kyryr r’k-ir-^-r-k-r-K-k'k'x-h-k-x 

MODES2 PROVIDES PRESS AKD DER AT 'LOV.ER' CI.'LY 
INPUTS 

VAL - THE INCREMENTAL HORIZONTAL WAVENUM U? TO 819 
DELZ- DEPTH INCREMENT EM OF DEPTH, t."JK3 CE MODES, 
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7ILE; tL-p FORTK.*:: R. 1 



C OF INCREMENTS / MODE 

C :i - NUMBER OF MODES 

C R - HORIZONTAL WAVENUMBER (COMPLEX) AT EACH DEPTH IN-C 

C DUP - DEPTH OF UPPER OF TXER AND RXER 

C DLW - DEPTH OF LOV.'ER OF TXER RXER 

C IU1/IU2 INCREMENT N'UMBER FOR INC ABOVE AND BELOW ’ Upp 

C IL1/IL2 INCREMENT NUMBER FOR INC ABOVE ANT BELOW 'LOW 

C INC - NUMBER OF INCREMENTS OF DEPTH 

C OUTPUTS 

C U1 - 'PRESSURE' AT UPPER 

C DUl - 'DERIVATIVE AT UPPER 

C VI - 'PRESSURE' AT LOWER 

C DVl - 'DERIVATIVE' AT LOWER 

C 

-V TV y; -k -k -r -k if ie -k "k -r ie -k rk -k -k -ir -k -r •ir if -k -K if if -k ir -x Tf -k -K ^ re ir Tf yc ^ -r y- -X yr -r -K ->r •k ^ -W -k -X -x -r 

c 

SUBROUTINE MODES2 { VAL , DELZ , INCTOT, K, V,DV, ILl, IL2.DLW) 
COMPLEX K ( 25000 ) , P ( 2 ) , DP ( 2 ) , FI , A1 , 5 1 , F2 , A2 , E2 , 

1 F3 , A3 , B3 , F4 , A , B , PD I FF , QD I FF , V , DV , VAL 
B = .00001 

A = 0.0 



DO 28 


L 


= 1,INCT0T^1,2 




LI 


= 


IMCTOT - L - 1 




FI 


= 


-B (K(Ll)--2 


- VAL) 


A1 


= 


A + DELZ-Fl 




B1 


= 


B + DELZ-A 




F2 


= 


-Bl*(K(Ll-l)--2 


- VAL) 


A2 


= 


A + DELZ’^F2 




B2 


= 


B + DELZ^^Al 




F3 




-B2*(R(Ll-l)*-2 


- VAL) 


A3 


= 


A + 2*DELZ*F3 




B3 


= 


B + 2-DELZ-A2 




FT 


= 


-B3^(K(Ll-2)-^2 


- VAL) 



B = B + 2’-DELZ-(A * 2-Al - 2-A2 + A3) '6 
A = A + 2^-DELZ’(Fi + 2^-F2 -*■ 2-^F3 F4)/6 

C 

IF( (Ll-2) ,EQ. ILl) THEM 
P(l) = B 
DP(1)= -A 
END IF 

IF((Ll-2) .EQ. IL2) THEN 
P(2) = B 
DP(2) = -A 
ENDIF 
C 

28 CONTINUE 

PD IFF = P(2) - P(l) 

QDIFF = DP(2) - D?(l) 

DDIFF = DLW - (ILl)-DELZ 
V = P(l) + PDIFF-DDIFF/(DELZ-2 ) 

DV = DP(1) + QDIFF*DDIFF/(DELZ-2 ) 

C 

21 RETURN 
END 
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